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Abstract 


A  method  of  solution  for  the  equation  of  radiative 
transfe  for  a  spherical,  grey  atmosphere,  steady  state 
plasma  in  radiative  equilibrium  is  developed.  The  method  is 
called  the  ray-integration  technique  and  derives  from  the 
same  method  of  solution  done  in  cylindrical  geometry  by 
George  Nickel  of  Los  Alamos  National  Laboratory.  The  total 
and  net  radiation  flux,  source  function,  opacity  and  moments 
of  intensity  are  calculated  as  a  function  of  radius.  The 
solution  to  the  zero' th  and  first  moment  equations  are  also 
provided . 

The  conditions  of  the  grey  atmosphere  problem  and  the 
general  nature  of  radiation  transport  in  spherical  geometry 
are  developed  and  the  ray  integration  technique,  as  applied 
to  the  problem,  is  presented. 

Numerical  results  are  calculated  for  a  variety  of  radial 
mesh  and  opacities.  These  results  are  provided  in  graphic 
form  and  compared  against  theoretically  predicted  behavior. 

The  ray-integration  technique  developed  in  this  thesis 
produces  numerical  results  which  are  in  reasonable  agreement 
with  theoretical  predicted  results. 


THE  RAY-INTEGRATION  TECHNIQUE 
IN  SPHERICAL  GEOMETRY 


I.  Introduction 


Background 

Motivation.  In  any  study  of  radiative  energy  transport, 
the  formulation  and  solution  of  a  fundamental  conservation 
equation  is  of  paramount  importance.  This  equation  is  known 
as  the  equation  of  radiative  transfer,  or,  more  commonly,  the 
transfer  equation.  The  transfer  equation  predicts  the  char¬ 
acter  of  the  radiation  field  of  a  plasma  by  solving  for  the 
intensity  at  a  point  within  the  field  as  function  of  time  and 
position.  This  thesis  prescribes  a  method  of  solution  to  the 
transfer  equation  in  a  one-dimensional,  steady-state,  fre¬ 
quency  independent  medium  in  one-dimensional  spherical  geom¬ 
etry. 

The  method  used  is  called  the  ray-integration  technique 
and  is  based  on  work  done  by  George  Nickel  of  Los  Alamos  Na¬ 
tional  Laboratory. 

TRAI LMASTER  Project .  The  ray-integration  technique  is 
derived  from  work  done  by  Nickel  in  support  of  the  U.  S.  Air 
Force's  TRAILMASTER  project.  The  goal  of  this  project  is  to 
create  a  high  energy  source  of  soft  X-rays  from  the  rapid 
collapse  of  a  cylindrical  foil.  In  support  of  this  goal,  Dr. 
Nickel  created  a  program  called  LMILNE  which  executed  the  ray- 
integration  technique  in  cylindrical  geometry.  The  program  he 


created  incorporated  some  techniques  which  could  be  expanded, 
with  modification,  to  different  geometries.  It  was  determined 
that  so  modifying  the  technique  into  spherical  geometry  would 
prove  useful. 

General  Applicability.  A  solution  for  the  transfer  equa¬ 
tion  is  frequently  needed  in  many  branches  of  physics.  The 
preponderance  of  dense  spherical  plasmas  with  radiation  fields 
modeled  in  research  provides  the  motivation  to  find  methods  of 
solution  which  are  both  efficient  and  accurate.  The  ray-in¬ 
tegration  method  has  proved  to  be  both. 

Analysis  of  Stellar  Atmospheres .  Applied  to  spherical 
geometry,  the  ray-integration  method  of  solving  the  equation 
of  radiative  transfer  could  be  very  useful,  especially  to  the 
astrophysics  community  as  a  "first  order  estimate"  device  for 
modeling  extended  stellar  atmospheres.  The  assumptions  made 
about  the  nature  of  the  plasma  system  solved  for  in  this  in¬ 
vestigation  is  completely  analagous  to  the  radiative  equilib¬ 
rium,  grey  atoinosphere  model  often  used  in  text  books  to  illus¬ 
trate  the  basic  nature  of  radiation  transport.  The  research 
for  this  thesis  depended  extensively  on  this  analogy.  Although 
normally  referring  to  a  stellar  body,  the  terms  atmosphere  and 
interior  are  used  throughout  this  paper  to  refer  to  the  regions 
of  energy  transport  and  source  region  respectively.  This  in  no 
way  alters  the  applicability  of  the  assumptions  and  solutions 
to  spherical  plasmas  in  general. 


Objectives  and  Scope 

The  goal  of  this  research  has  been  to  produce  a  physi¬ 
cally  realistic  solution  to  the  transfer  equation  which: 

--  produces  numerically  accurate  results  from  which  the 
character  of  a  plasma  can  be  approximated 

--  can  be  used  to  engender  an  intuitive  understanding  of 
the  physical  nature  of  stellar  atmospheres. 

A  second  objective  of  this  research  is  to  verify  the  ray- 
integration  technique  as  plausible  by: 

--  implementing  a  computer  program  which  utilizes  this 
technique 

--  producing  graphs  of  the  results  using  various  combina¬ 
tions  of  radii,  opacity,  power  and  iteration  techniques 

--  comparing  the  results  of  this  program  against  theoret¬ 
ically  predicted  results. 

Assumpt ions 

To  effectively  develop  a  method  of  solution  for  the  be¬ 
havior  of  a  system,  its  physical  character  must  be  established 
That  is,  assumptions  about  the  size,  shape,  and  material  prop¬ 
erties  of  the  system,  in  this  case,  the  plasma,  must  be  known 
beforehand.  This  does  limit  the  applicability  of  a  technique 
to  a  specific  physical  model,  but  often  provides  sufficient  in 
sight  to  generalize  the  method  to  other  systems. 

The  particular  assumptions  made  are  that  the  ray-integra¬ 


tion  technique,  as  developed  in  this  thesis,  is  applicable 
specifically  to  a  system  that  is 


--  frequency  independent  (called  the  grey  atmosphere 
problem) 

--in  radiative  equilibrium 

--  time  invariant,  i.e.  steady  state. 

Method  of  Solution 

The  ray-integration  technique  is  implemented  using  a  com¬ 
puter  program  written  in  the  FORTRAN  language.  The  program 
was  executed  on  both  a  CRAY  located  at  Los  Alamos  National 
Laboratory  and,  with  minor  changes,  on  the  author's  APPLE  11+ 
compatible  computer.  The  graphs  shown  in  this  report  were 
done  exclusively  on  the  CRAY. 

The  inputs  to  the  program  were  values  for  the  opacity  and 
source  function  of  the  radiation,  as  well  as  the  number  of  ra¬ 
dii  and  radii  spacing  which  consequently  made  up  the  radii 
mesh . 

The  output  of  the  program  is  the  net  radiation  flux,  the 
total  radiation  flux,  and  the  ratio  of  moments  of  intensity 
known  as  the  Eddington  Factors.  Additional  output  consists  of 
solutions  to  the  moments  of  the  transfer  equation.  All  input 
and  output  are  in  MKS  units. 

Plan  o f  Presentation 

This  report  is  presented  in  the  following  manner.  The 
first  chapter  is  a  review  of  the  basic  theories  of  radiation 
transport  and  of  the  specific  nature  of  the  grey  atmosphere 
problem  in  radiative  equilibrium.  Chapter  II  presents  an 


overview  and  develops  the  ray-integration  technique.  The  next 
chapter  provides  a  development  of  the  numerical  methods.  This 
is  followed  by  an  outline  of  the  program  developed  to  imple¬ 
ment  the  ray-integration  technique.  Next  presented  is  a  sum¬ 
mary  of  some  results  generated  by  the  computer  code,  along  with 
a  discussion  of  these  results.  Finally,  conclusions  and  recom¬ 
mendations  for  further  study  are  put  forth  in  the  last  chapter. 


II .  Theory 


To  understand  the  technique  for  solving  the  equation  of 
radiative  transfer,  it  is  first  necessary  to  review  some  of 
the  basic  physical  theories  of  the  transfer  of  energy  and  de¬ 
velop  the  transfer  equation  for  analysis.  This  will  be  done 
in  this  section  beginning  with  the  characteristics  of  the  ra¬ 
diation  field  and  progressing  through  a  general  description 
of  the  grey  atmosphere.  The  development  of  the  following  con¬ 
cepts  is  fairly  standard  and  can  be  found  in  many  textbooks 
such  as  references  2,  3,  6,  7,  and  9. 

Description  of  the  Radiation  Field 

The  region  within  and  surrounding  the  body  of  a  star  or 
spherical  plasma  mass  is  a  radiation  field  which  contains  pho¬ 
tons  of  various  frequencies  traveling  in  all  directions.  The 
interaction  of  the  photons  with  the  mass  of  the  body  and  with 
each  other  results  in  energy  transfer.  A  description  of  the 
intensity  variations,  then,  describe  the  nature  of  the  radia¬ 
tion  field  from  which  other  material  properties  of  the  mass 
may  be  derived. 

Specific  Intensity.  The  radiation  field  produced  by  the 
emission  and  absorption  of  photons  brings  about  a  transport  of 
radiant  energy,  dE,  in  a  specified  energy  interval  (v,  v  +  dv). 
Picture  this  energy  transfer  as  a  beam  of  radiation  flowing 
from  a  surface  element  area  dA,  into  a  solid  angle  dfi ,  in  a 
direction^  that  makes  an  angle  (<p)  with  the  normal  to  the 
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surface,  A,  (see  Figure  1). 


Figure  1 

A  beam  of  radiation,  originating  from  surface 
element  dA,  into  a  solid  angle  dft, 
travelling  in  direction  ft  of  length  s. 

The  transfer  of  energy  along  a  distance  r  and  through 
time  interval  dt  is  expressed  by  the  equation 

dE(n,j3,vit)  =  I(n,£,v,t)  dft  dt  dv  cos(cp)  ds  (1) 

where  I(n,£,v,t)  is  defined  as  the  specific  intensity  (or, 

more  simply,  intensity).  Suppressing  the  time  and  spatial 

dependence,  intensity  =  I  .  The  units  of  intensity  are  watts 
-2  -l  -l 

m  sr  hz  It  is  clear  that,  by  knowing  the  specific  in¬ 

tensity  at  points  progressing  through  a  region,  the  energy 
present,  and  many  quantities  derivable  from  the  energy,  can 
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Moments  of  Intensity 

The  Mean  Intensity.  In  addition  to  intensity,  there 
are  other  quantities  used  frequently  to  describe  the  behavior 
of  the  radiation  field.  These  are  the  angular  averages,  or 
the  angular  moments,  of  intensity.  There  are  three  moments 
which  are  used  in  spherical  geometry.  The  first  is  called  the 
zero' th  order  moment  or  the  mean  intensity,  J,  and  is  simply 
the  average  of  the  intensity  over  all  angles  (thus,  it  is  of¬ 
ten  also  called  the  "all-angle  intensity").  The  mean  inten¬ 
sity  is,  therefore,  given  by 


J  -  t~  fn  I  dn 
v  4  it  J  ft  v 


(2) 


The  units  of  J  are  the  same  as  the  units  for  I  . 

v  v 


It  is  common  practice  to  transform  the  integral  over 
solid  angles  to  an  integral  over  the  quantity  y  =  cos((p). 
This  is  done  by  substituting  the  differential  dft  by 
sin(cp)  d&  dtp  =  -  d$  dy.  If,  as  is  usually  the  case,  there 
is  symmetry  about  the  azimuthal  direction,  then  the  integral 
becomes 


Jx.  =l/_,Iv(u)  du 


(3) 


Expressing  the  angular  moments  in  terms  of  integrals  in  m 
is  known  as  Eddington  notation. 


The  Flux.  The  first  order  moment  of  intensity  is 
known  as  the  flux,  £  .  Physically,  this  flux  represents  a 
vector  quantity  such  that  (5  •  dA)  gives  the  radiation  en¬ 

ergy  flow  across  a  surface  area  dA  per  unit  time,  per  unit 
frequency  interval.  Analytically,  this  is  given  by 

Cv  =  /flIv(£i)  cos(cp)  dft  (A) 

In  Eddington  notation,  this  integral  becomes 

=  2u  J  ^  (p)  \i  d y  (5) 

-  2  -  1 

The  units  of  total  radiation  flux  are  watts  m  hz  .  Addi¬ 
tionally,  a  quantity  called  the  Eddington  flux  is  defined  as 

Hv  =  ?/',  Vu>  u  du  (6) 

It  can  be  seen  that  the  Eddington  Flux  is  an  average  flux 
per  steradian.  Its  units  are  watts  m  2str  1  hz  1  . 

One  other  commonly  seen  version  of  the  radiation  flux  is 
called  the  astrophysical  or  total  flux,  F.  This  is  simply 
the  flux  emitted  over  all  directions  and  is 


F  (r)  =  Anr2  H 
v  v 


(7) 


where  r  is  the  distance  from  the  center  of  the  sphere  to  the 
location  at  which  the  flux  is  being  evaluated. 


The  Second  Moment .  The  second  moment  of  intensity, 

K,  is  related  to  the  radiation  pressure  (3:68).  This  quantity 

2 

is  simply  constructed  by  multiplying  the  intensity  by  cos  (<p) 
and  integrating  over  all  solid  angles.  Using  Eddington  nota¬ 
tion  , 

K  =  j  f1  I  (y)  u2  dy  (8) 

v  1  J  -i  v 

The  units  for  K  are  the  same  as  for  I  ,  J  ,  and  H  . 

v  v  v 

Summary.  The  traditional  radiation  quantities 
I,  J,  H,  and  K,  described  in  the  last  sections  are  generally 
sufficient  to  develop  an  approximate  analysis  of  the  radia¬ 
tion  field.  For  emphasis  and  easy  referral,  they  are  summa¬ 
rized  in  Table  1.  The  equations  are  given  in  Eddington  no¬ 
tation  because  this  the  form  will  be  used  throughout  the  re¬ 
mainder  of  this  report. 


Table  I 

Summary  of  Intensity  and  Moments  of  Intensity 
Nomenclature  Symbol  Equation 


Intensity 

Mean  Intensity 
(Zero  Moment) 


Eddington  Flux 
(First  Moment) 


I  dE  =  I  dft  dt  dv  cos(cp)  ds 

v  v 

J  J  =  i  /  I  (y)  dy 

V  V  Z  J  -  1  V 


H  H  =  i  /  I  (y)  y  dy 

v  v  2  J  _  i  v 


kv  *  \  C  V">  du 


Second  Order 
Moment 


The  Interaction  of  Radiation  with  Matter 

Having  thus  defined  the  quantities  which  approximate  the 
radiation  field,  we  now  consider  the  processes  which  produce 
and  modify  them. 

Interaction  Processes .  Photons  are  catagorized  into  spe¬ 
cific  energy  groups  of  given  frequency  intervals.  The  pho¬ 
tons,  as  they  pass  through  the  matter  in  the  radiation  field, 
constantly  interact  with  that  matter  and  with  each  other. 

They  may  be  scattered  out  of  an  energy  interval  or  absorbed 
by  atoms,  ions  or  molecules.  They  may  be  emitted  into  energy 
groups  by  spontaneous  decay,  or  simply  scattered  into  an  en¬ 
ergy  group  as  a  result  of  being  scattered  out  of  another.  The 
result  of  these  interactions  is  to  modify  the  intensity  of  the 
radiation  field.  Detailed  descriptions  of  these  processes  are 
not  necessary  for  the  comprehensive  development  of  this  thesis. 
It  is  sufficient  to  separate  these  interaction  processes  into 
two  groups,  removal  mechanisms  and  addition  mechanisms,  and 
consider  only  the  microscopic  results  of  the  processes  in  terms 
of  opacity  and  emissivity. 

Opacity .  The  probability  that  a  photon  will  be  either 
scattered  or  absorbed  in  a  prescribed  path  length  is  given  by 
the  extinction  coefficient,  or  opacity,  xv*  The  opacity,  in 
general,  is  defined  in  a  similar  manner  as  the  intensity.  That 
is,  the  amount  of  material  removed  can  be  represented  by  an  en¬ 
ergy  flow,  such  that  a  volume  of  material  of  cross  section  dA, 
and  length  da,  removes  from  a  beam  with  specific  intensity  I  , 
incident  normal  to  dA  and  propagating  into  a  solid  angle  dft,  an 


amount  of  energy 


6E(n,sj,v,t)  +  I(n,£,v,t)  dA  da  df2  dt  (9) 


within  a  frequency  interval  dv,  in  a  time  dt  (5:23). 

There  are  different  opacities  for  scattering  and  for  ab¬ 
sorption,  and  even  different  opacities  for  the  categories  of 
scattering  (elastic,  coherent,  etc.)  and  for  the  types  of  ab¬ 
sorption  that  occur  ( photonization ,  free-free  absorption, 
etc.)  However,  since  these  processes  are  statistically  inde¬ 
pendent,  the  opacities  add  linearly  such  that 


(10) 


where  y  s  is  the  scattering  opacity,  y  a  is  the  absorption 
T 

opacity  and  y^  is  the  total  materials  opacity.  For  conven¬ 


ience,  in  this  thesis,  no  distinction  will  be  made  between 
the  particular  type  of  removal  mechanism  and  the  term  opacity 
will  be  used  to  refer  to  the  total  material  opacity.  Since 
opacity  is  the  probability  that  a  photon  will  interact  over 
a  given  distance,  the  inverse  of  the  opacity  is  simply  the 
mean  distance  that  a  photon  will  go  before  an  interaction  oc¬ 
curs,  i.e.,  the  inverse  of  the  opacity  is  the  mean  free  path. 
The  units  of  opacity  are  length  1 ;  or  m  1 . 

Emissivity .  The  population  of  photons  in  a  frequency 


interval  may  be  altered  not  only  by  removal  from  the  interval, 
but  also  by  addition  to  the  interval  from  other  frequency 


m 

m 


groups.  This  may  be  done  either  by  the  photons  being  scattered 
into  the  group  as  they  were  scattered  out  of  others,  or  they 
may  be  created  by  some  emission  process. 

Following  the  process  used  to  define  the  extinction  coef¬ 
ficient,  we  define  the  emission  coefficient,  or  emissivity, 
n  ,  such  that  the  amount  of  energy  released  from  an  element  of 
material  of  cross  section  dA  and  length  da  into  a  solid  angle 
dft,  within  a  frequency  band  dv ,  in  direction  £  in  a  time  inter¬ 
val  dt,  as 


3E(n,s^\>,t)  =  ny(£)S,\))t)  dA  da  dft  dv  dt 


(ID 
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The  units  of  emissivity  are  watts  m  sr  hz  (7:25).  Again, 
the  above  derived  quantity  is  the  total  emissivity  and  is  a 
linear  sum  of  emissivities  of  the  separate  emission  processes. 

Optical  Depth  and  the  Source  Function 


Having  defined  the  concept  of  opacity  and  emissivity,  two 
very  important  concepts  are  presented  which  derive  from  the 
above.  These  are  the  optical  depth  and  the  source  function. 

The  Optical  Depth .  Given  the  opacity  along  a  path  length 
ds,  the  differential  optical  thickness  is  defined  such  that 


d  t  =  y  ds 
v 


(12) 


Integrating  over  a  distance  along  the  line  of  sight  yields 


t  (s)  =  /  2  x  (s') 


(13) 
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where  is  now  the  optical  depth.  Since  ds  is  in  units  of 
distance  and  y  is  in  units  of  inverse  distance,  t  is  dimen- 
sionless.  The  limits  si  and  s2  are  defined  such  that  the  in¬ 
tegral  is  always  positive,  that  is  the  optical  depth  increases 
as  r  ->  0  and  is  0  at  the  surface  of  the  sphere.  The  optical 
depth  is  the  number  of  mean  free  paths  along  the  line  of  sight. 
A  region  is  said  to  be  "optically  thick"  if  x^  >>  1,  and  "op¬ 
tically  thin"  if  x  <<  1. 

v 

The  Source  Function .  As  stated  in  section  B,  the  radia¬ 
tion  field  is  continually  being  perturbed  by  the  removal  and 
addition  of  photons.  The  ratio  of  the  amount  added  to  the 
amount  removed  is  simply  the  ratio  of  the  emissivity  to  the 
opacity  and  is  called  the  source  function,  S.  The  source  func¬ 
tion  is,  then,  given  by 


S  = 
v 


(14) 


where  is  the  total  opacity  and  nv  is  the  total  emissivity. 
Although  simply  defined,  the  source  function  can  actually  be 
quite  involved  to  calculate,  depending  on  the  number  of  proc¬ 
esses  occurring  which  produce  the  opacity  and  emissivity  of 
the  material.  An  excellent  treatment  of  some  of  the  more  com¬ 
mon  variations  of  the  source  function  can  be  found  in  Refer¬ 
ence  7,  Chapter  2  and  Reference  2. 

The  Transfer  Equation 

The  material  properties,  opacity  and  emissivity,  discussed 
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in  Section  C,  modify  the  intensity  of  a  stellar  atmosphere. 
The  change  produced  can  be  expressed  in  words  as 


[change  in  I  per  path  length  ds]  =  [source  -  sinks]  (15) 


where  the  term  "sources"  means  change  in  intensity  per  path 
length  ds  due  to  the  emission  processes  occurring,  and  "sinks" 
refers  to  the  changes  occurring  which  are  due  to  absorption. 
Thus,  the  above  can  be  rewritten  in  terms  of  opacity  and  emis- 
sivity  as 


dl  _  T 

ds  nv  Vv 


(16) 


where  xv  is  tlie  total  materials  opacity  and  is  the  total 
materials  emissivity.  Expressing  the  left  hand  side  explic¬ 
itly  in  terms  of  position  and  time  yields 


dl  _  1  31  .  n  nT 
v  -  —  v  +  ft  'VI 


(17) 


Coupling  equations  (16)  and  (17), 
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which  is  the  transfer  equation  in  differential  form.  Using 
(14)  and  assuming  a  steady-state  radiation  field, 


ft. VI  =  x  (S  -  I) 

—  V  V  V 


(19) 


. 


mm®. 


For  a  one-dimensional,  spherically  symmetric  radiation 
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field, 


/  n  31  sm((p)  31  _  ( c  T  \ 

cos(o)  -t—  -  — — *—  -r—  -  y  (.S  -  I  ; 
Y  3r  r  3cp  Av  v  v 


3cp 


(20) 


In  terms  of  Eddington  notation,  this  is 


41  +  —  4- ■■■  41  =  y  (S.  -  I  .) 


3r 
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(21) 


Yet  another  form  of  the  equation  (19)  is  produced  using  the 
definition  of  the  source  function  and  the  relationship 


dx  =  v  dS ,  such  that 
v  Av 


dl 


v  =  S  -  I 

cTT"  v  V 


(22) 


Multiplying  both  sides  of  the  above  by  ,  and  integrating 


over  dx,  from  x  =  0  to  x  =  x  .  The  result  is 
v  v  v 


-x 


T  (X  ’ _  X  ) 

I  (x  )  «  I  (o)  e  v  +  J  VS  e  v  y  dx  1 

V  V  V  V  V 


(23) 


A  number  of  variations  of  the  transfer  equation  have  been 
presented  in  this  section.  To  review,  and  for  the  convenience 
of  the  reader,  the  following  table  lists  the  version  of  the 
transfer  equation  and  its  solution  which  will  appear  explic¬ 
itly  in  a  later  section. 
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Table  II 

The  Transfer  Equation  and  its  Solution  for  a  One-Dimen¬ 
sional,  Spherically-Symmetric  Radiation  Field 


Name 


Equation 


Transfer  eq. 
in  spherical 
coordinates 


31 

> 

3y 


=  x  ( s 

v 


I 

V 


) 


Solution  to  the 
transfer  eq. 


= 


Iv(o)e  Tv 


x' 

v 


The  Moments  Equations 

As  the  transfer  equation  represents  the  change  in  inten¬ 
sity,  so  the  moments  equations  represent  changes  in  the  mo¬ 
ments  of  intensity.  These  equations  will  now  be  presented. 

The  Zero ' th  Moment  Equation .  The  zero' th  moment  equation 
is  formed  from  equation  (21).  Integrating  this  equation  over 
all  solid  angles  or,  in  Eddington  notation  integrating  over  p 
from  -1  to  1 ,  yields 


d  u  + 


1  r1  (1-u2) 

2  - 1  r 
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-  I  )  d u  (24) 
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Using  the  relation 
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Equation  (24)  becomes 
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|f[  i  L,  M  \  dy  1  +  r  ["  \  L,  u  *v  dul  =  “T” -L^  +TL[(.I1Iv  dul  (26) 


Substituting  in  for  the  boxed  quantities,  those  moments  they 
equate  to  (see  Table  I,  page  10)  the  zero'th  moment  equation 
is  revealed  to  be 


^  +  -  H  =-  =  x  (S  -  J  ) 

3  r  rvr2  9  r  v  v 


(27) 


The  First  Moment  Equation.  Proceeding  in  the  same  manner 
as  above,  equation  (21)  is  this  time  multiplied  through  by  y 
and  integrated  over  all  y.  Again  using  relation  (25)  produces 


I?  [ l(!,  V2  TV  du  1  +  7  ['  dy  +  hj  »2  K  dM  1  =  [/l1Iv  dw  1 


(28) 


Again,  substituting  the  appropriate  moments  from  Table  I 
yields  the  first  moment  equation. 


9K  +  3  K  -  J 
3r  r 


-  x  H 

v 


(29) 


Summary .  For  emphasis,  the  following  table  lists  the 
significant  equations  of  this  section. 
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Table  III 


The  Moment  Equations 


Name 

Equation 

Zero ' th 

Moment  Eq . 

,  / r 2H  . 

1  vl  _  Y 

r2  9r 

<sv  -  V 

(27) 

First 

9Kv  .  3Kv_Jv 

=  -XvHv 

(29) 

Moment  Eq. 

3r  +  r 

Assumptions  and  Approximations 

Vital  to  any  scientific  work  are  the  assumptions  and  ap¬ 
proximations  that  are  made.  Already  certain  approximations 
have  been  presented  in  this  thesis.  In  Section  C,  the  approx¬ 
imate  regions  traditionally  considered  optically  thick  or  thin 
were  outlined.  In  Section  B,  causally  different  inter-atomic 
processes  were  grouped  under  the  same  headings  as  either  re¬ 
moval  or  addition  mechanisms.  The  transfer  equation  in  Sec¬ 
tion  D  was  specialized,  based  on  a  steady  state  approximation. 
All  these  are  valid  in  general  and  essential  to  the  develop¬ 
ment  of  this  thesis.  In  this  section,  the  last  approxima¬ 
tions  and  assumptions  necessary  to  be  stated  are  presented. 

The  Eddington  Approximation .  The  material  properties, 
opacity  and  emissivity  modify  the  energy  flow  and,  thus,  the 
intensity  and  its  moments.  In  particular,  within  certain 
limits  of  optical  thickness,  the  moments  take  on  approximate 
values.  This  section  discusses  one  very  frequently  used  ap- 
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proximation,  the  Eddington  approximation,  and  variations  of  it. 


•V 


The  Eddington  Factor . 

The  Diffusion  Limit .  The  Eddington  approxima¬ 
tion  is  manifested  in  the  variable  Eddington  factor,  f^,  which 
is  defined  as  the  ratio  of  the  second  moment  of  intensity  to 
the  zero' th,  or,  more  simply,  f^  =  K/J.  The  approximation  is 
that,  at  great  optical  depths,  the  radiation  field  is  nearly 
isotropic.  The  area  where  this  is  true  is  called  the  diffusion 
limit.  In  this  limit,  the  intensity  is  subsequently  independ¬ 
ent  of  angle  and  can  be  removed  from  the  moments  integrals  of 
Table  II,  page  17.  Doing  so  reduces  those  equations  to 

J  =  I  i  f1  dp  =  I  (30) 

v  \i2Lj  v 

and 

i  i  1 

K  =  I  i  /  y2  d u  =  -^  (3D 

Av  V  2  J  _  1 M  M  3 

Thus,  in  the  diffusion  limit,  =  1/3.  This  is  physically 
valid  and  can  be  derived  independently  as  part  of  the  diffu¬ 
sion  approximation  (7:51). 

The  Streaming  Limi t .  While  the  Eddington  ap¬ 
proximation  certainly  holds  true  in  the  optically  thick  re¬ 
gions  of  the  plasma,  what  then  can  be  said  about  its  accuracy 
in  the  optically  thin,  less  dense,  outer  regions?  In  this 
case,  the  Eddington  Approximation  must  be  augmented. 

Consider  that,  in  the  outer  region,  the  photons  are  no 
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longer  colliding  as  frequently,  since  their  population  is,  by 
definition,  sparse.  Thus,  the  mean  free  path  is  comparatively 
long  and  the  photons  can  be  thought  of  as  essentially  travel¬ 
ing  in  uninterrupted  paths  or,  streaming.  In  this  streaming 
limit,  all  photons  are  going  in  the  same  "direction":  out. 
Therefore,  all  angular  moments  must  sum  to  the  same  value  and 
any  ratio,  thereof,  must  be  one,  thus  the  Eddington  factor, 
f ^ ,  is  one  in  the  streaming  limit.  It  is,  therefore,  more  ac¬ 
curate  to  talk  about  not  an  Eddington  factor,  but  a  variable 
Eddington  factor:  one  which  varies  from  1/3  at  deep  optical 
depths  within  the  system  to  1  at  the  outer  boundry  of  the  sys¬ 
tem. 

The  "Second"  Eddington  Factor .  Although  most  liter¬ 
ature  refers  only  to  the  ratio  f^  as  the  Eddington  factor, 
there  is  also  another  comparative  ratio  which  is  frequently 
used.  This  is  the  ratio  of  the  Eddington  flux  to  the  mean  in¬ 
tensity,  or  H/J.  This  "second"  Eddington  factor  will  be  re¬ 

ferred  to  in  this  thesis  as  fH.  For  convenience,  both  ratios 
will  be  termed  the  Eddington  factors  and  will  be  distinguished 
by  their  subscripts.  As  with  the  first  Eddington  factor,  the 
ratio  H/J  also  has  limiting  values. 

The  Flux  a£  £  _£>  0.  In  the  central  regions  of 

the  sphere,  it  is  clear  to  see  that  fjj  must  be  zero.  This  is 

due  to  the  fact  that  as  r  ->  0,  the  radiation  field  is  iso¬ 
tropic  and,  thus,  the  total  radiation  flux  must  be  zero,  since 
there  is  no  net  energy  flow.  Analytically,  this  can  be  seen 
by  returning  to  Table  II  and  noting  that  for  an  isotropic 
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field  I  is,  by  definition,  independent  of  and  can  be  removed 
from  the  integral.  What  remains  then  is 

H  =  I  i  J  *  y  dp  =  0  (32) 

v  v  ^  - 1 

Thus,  if  H  ->  0,  as  r  ->  0,  so  .  to  does  fjj. 

The  Streaming  Limi t .  In  the  streaming  limit 
fH  ->  1.  This  is  so  for  the  same  reasons  given  previously 
for  the  streaming  limit  approximation  of  f^. 

Radiative  Equilibrium.  A  stellar  atmosphere  by  defini¬ 
tion  is  that  portion  of  the  star  that  transmits  the  radia¬ 
tion  which  is  created  deep  within  the  interior.  In  the  at¬ 
mosphere  region,  these  energy  transfer  processes  take  on  two 
forms  in  general:  radiative  and  hydrodynamic.  For  the  pur¬ 
poses  of  this  thesis,  it  is  assumed  that  all  transfer  proces¬ 
ses  are  radiative.  In  other  words,  it  is  assumed  that  the 
plasma  system  is  in  radiative  equilibrium;  that  is,  all  en¬ 
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ergy  is  transferred  through  radiative  processes,  thus  all  en¬ 
ergy  absorbed  in  a  volume  element  must  equal  that  radiated  in 
that  element.  If  over  all  frequency  intervals,  the  energy  re¬ 
moved  from  the  beam  is 


Jo”  d v  <j>  dn  xv  Iv  =  4°°  xv  dv 


(33) 


And  the  energy  emitted  is 


J  dv  <£  d  ft  n  =  4  n  J  x  S  dv 

Jo  '  y  •'®AVV 


(34) 
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or,  J  =  S  ,  (7:49).  Note  that  the  radiative  equilibrium 
v  v 

assumption  states  that  what  is  created  within  a  volume  is  de¬ 
stroyed  in  this  volume  implies  that  the  total  radiation  flux 
within  a  volume  must  be  constant,  or 
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or 


r2  H  =  constant:  4nr2H  =  F  =  constant  (36) 

v  v  v 

The  Grey  Atmosphere .  A  fundamental  example  for  the 
analysis  of  stellar  atmospheres,  and  the  specific  case  of 
study  of  this  thesis,  is  the  investigation  of  the  grey  at¬ 
mosphere.  The  grey  atmosphere  problem  is  appropriate  for 
the  study  of  solutions  to  the  transfer  equation  because  the 
assumptions  which  are  entailed  are  such  that  the  solution 
of  the  equation  becomes  independent  of  the  nature  of  the  ma¬ 
terial  itself.  In  addition,  the  problem  introduces  the 
aspect  of  radiative  equilibrium,  which  may  then  be  expanded 


to  more  realistic  cases  of  non-grey,  non-static  problems. 
The  crux  of  the  grey  atmosphere  problem  is  that  the 


opacity  of  the  material  is  independent  of  frequency,  thus 

=  x  (a  constant).  Further,  integrating  all  quantities  over 
frequency  such  that 


(37) 


where  x  =  I  ,  J  ,  S  ,  etc.,  the  transfer  equation  takes  on 

v  v  v  v 

the  form 


u  il  +  Il.~ M 2) -II  =  y  (S  -  J)  (38) 
3  r  r  3p  * 

Summary .  Before  finally  moving  on  to  the  specifics  of 
the  ray-integration  technique,  a  summary  of  the  assumptions 
and  approximations  of  this  section  is  given  below. 


Table  IV 

Radiative  Equilibrium  and  Grey  Atmosphere  Assumptions 


Assumption 


Statement 


Result 


Radiative  J  =  S  4  irrzH  =  constant 

Equilibrium 

Grey  X  =  I»  J>  H>  K>  s  are 

Atmosphere  frequency 

independent 


III.  The  Ray-Integration  Technique 


Having  developed  the  basic  theory  of  radiation  transport 
attention  is  now  turned  to  the  actual  solution  of  the  grey  at 
mosphere  problem  via  the  ray-integration  technique. 

Statement  of  the  Problem 

To  reiterate,  the  purpose  of  this  thesis  is  to  explore  a 
method  of  solution  to  the  spherical  coordinate  form  of  the 
equation  of  radiative  transfer,  specifically  for  the  grey  at¬ 
mosphere,  steady  state  problem.  This  method  of  solution  is 
called  the  ray-integration  technique.  As  the  name  implies, 
the  gist  of  this  method  is  the  solution  of  the  transfer  equa¬ 
tion  along  lines  of  intensity,  called  rays.  The  rays  are 
divided  into  small  length  segments,  ds,  and  the  transfer 
equation  solved  each  segment.  The  solutions  are  then  inte¬ 
grated  numerically  over  the  entire  ray.  A  step-by-step  de¬ 
velopment  of  this  method  is  described  in  the  next  section. 

Development  of  the  Method 

To  begin  an  explanation  of  this  method,  first,  imagine 
a  spherical  plasma  having  a  central  region  which  produces  co¬ 
pious  amounts  of  photons.  A  physical  analogy  of  this  would 
be  the  interior  of  a  star  where  subatomic  interactions  result 
in  large  amounts  of  radiation  output.  In  this  region,  radi¬ 
ation  field  is  very  dense  and  the  mean  free  paths  of  the  pho¬ 
tons  are  very  small.  In  other  words,  the  optical  depth  is 
very  large.  Here  the  diffusion  limit  is  valid.  Now  imagine 


that  this  interior,  very  dense  region  is  surrounded  by  an  ex 
tended  envelope  of  material  that  is  considerably  less  opaque 
to  photons,  i.e.,  a  region  where  the  photons  mean  free  path 
is  considerably  longer.  This  region  can  be  likened  to  the 
atmosphere  of  a  star.  It  is  through  this  atomosphere  region 
that  the  radiation  is  transported  from  the  interior  to,  and 
out  of,  the  surface  of  the  plasma  body. 

The  location  of  any  point  within  this  sphere  is  defined 
by  specifying  its  radial  distance ,  r,  from  the  center  and  the 
longitudinal  and  azimuthal  angles,  (a  and  8)  made  with  the 
axis'  of  the  sphere,  see  Figure  2. 


Figure  2 

Location  of  a  Point  on  a  Sphere 


Now,  subdivide  the  sphere  into  a  series  of  concentric  shells 
and  imagine  a  ray  of  photons  passing  through  the  sphere  tan¬ 
gent  to  one  of  the  shells.  The  position  of  the  ray  at  that 
tangent  point  is  now  further  defined  by  specifying  the  local 
or  directional,  longitudinal  and  azimuthal  angles  made  with 
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(dr)2  +  (  —  rdtp ) 2  »  ds2  (43) 

or 

d<£  =  -sin(  <p )  (44) 

ds  r 

Relating  these  relations  to  the  directional  derivative  used 
in  solving  the  transfer  equation,  we  have 
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(45) 


where  u  =  cos(<p).  The  above  is  just  the  relation  that  was 
simply  stated  in  the  last  chapter.  Recall  that  the  solution 
to  the  transfer  equation  is 


I  (x)  =  I  (o)  e-T  +  JT  S ( x )  e(T'"T)dx'  (23) 

which,  for  a  differential  length  of  ray,  ds,  and  using 

S2 

dT  =  x(s)ds>  x  =/slx(s')  ds',  is 
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The  object  of  the  ray-integration  technique  is  to  solve  the 
above  equation  for  each  length  of  ray,  ds,  between  two  radii 
and  to  add  the  solutions  for  each  increment  in  order  to  find 
the  solution  for  the  intensity  along  the  entire  ray;  in  other 
words,  to  numerically  "integrate"  the  intensity  over  the  en¬ 
tire  ray. 


Numerical  solutions  to  the  transfer  equation  via  the 
ray-integration  technique  were  implemented  explicitly  by  com¬ 
puter  code.  In  this  section,  the  numerical  methods,  the  logic 
behind  them  and  the  differencing  schemes  produced  are  detailed. 

Geometry  Dependent  Variables 

Radial  Mesh .  The  first  task  in  solving  a  geometry  depen¬ 
dent  problem  is  to  define  the  geometry  dependent  variables. 
First  of  these  are  those  variables  defining  the  physical  size 
of  the  plasma  and  the  number  and  spacing  of  radii  which  the 
sphere  is  divided  into;  in  other  words,  the  radial  mesh. 

The  size  of  the  sphere  is  defined  indirectly  by  inputting 
the  number  and  spacing  of  radii.  These  values  are  essentially 
arbitrary  and  can  be  scaled  up  to  the  size  of  a  stellar  body 
or  down  to  a  small  spherical  plasma.  The  form  of  the  radial 
mesh  is  also  arbitrary.  That  is,  an  evenly  spaced  mesh  can  be 
used  as  reasonably  as  a  logarithmic  scheme,  or  even  a  com¬ 
pletely  random  spacing.  An  example  of  how  a  five  radii, 
evenly  spaced  mesh  might  look  is  illustrated  by  Figure  A.  The 
mesh  is  indexed  such  that  jn  is  the  total  number  of  radii, 
with  the  first  radius  being  defined  as  zero.  That  is,  in  terms 
of  the  example,  jn  =  5,  r(l)  =  0.0,  r(2)  =  0.005  m,  etc., 
through  r(5)  =  0.02  m. 
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Figure  A 

Example  of  an  Evenly-Spaced  Radius  Mesh 


Angular  Mesh .  Refer  back  to  Figure  A.  Note  that  there 
is  an  angle  cp  made  with  the  radius  normal  for  each  ray  that 
passes  through  the  sphere.  If  many  rays  pass  through  the 
sphere,  all  at  different  cp  angles,  then  there  will  be  a  dis¬ 
tribution  of  angles  for  each  radius  so  that  each  radius  has 
an  angular  mesh  which  looks  like  Figure  5. 


Figure  5 

Angular  Mesh  of  an  Evenly-Spaced,  Five-Radius  Sphere 


We  see  that  each  radius  point  has  as  many  angles  associated 
with  it  as  rays  that  pass  tangent  to  it.  If  we  define  the 
rays  to  be  symmetric  about  p  =  0,  then  only  half  the  number 
of  angles  which  exist  need  to  be  generated,  since  the  angles 
on  one  side  of  the  centerline  are  just  the  negative  of  their 
mirror  mate. 

Generation  of  Angular  Weights.  Recall  that  the  quanti¬ 
ties  which  are  used  to  analyze  a  spherical  radiation  field  in¬ 
clude  the  angular  moments  of  intensity,  J,  H,  and  K.  Take  for 
example  the  moment  J  whose  equation  is 
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J  =  1  /  I (p )  dp 
7  p- 


(47) 


Where  p+  is  the  angle  at  the  beginning  of  the  ds  interval  and 
p-  is  the  angle  at  the  end.  This  integral  equation  must  be 
given  a  numerical  equivalent  for  solution  by  computer  code. 
This  was  accomplished  by  the  method  of  Gaussian  quadrature, 
such  that 


i  JU+Ku)  d v  *  l  I(u)  w. 


(48) 


where  w  are  the  numerical  weights  which  allow  the  summation 
to  approximate  the  integral  (see  Appendix  A  for  an  explana¬ 
tion  of  Gaussian  quadrature).  As  before,  I  is  the  intensity 
on  the  As  interval.  Note  that  intensity,  I,  must  be  given  a 
numerical  representation  as  well.  Thus  far,  the  intensity 


has  been  defined  in  terms  of  energy  flow,  but  no  particular 
polynomial  representation  has  been  presented.  There  are  actu¬ 
ally  a  number  of  ways  in  which  the  intensity  can  be  expanded, 
all  with  equally  varied  degrees  of  accuracy.  Many  forms  of 
polynomial  expansion  were  attempted  in  the  solution  of  this 
problem.  It  was  found,  however,  that  a  second  degree  Legendre 
polynomial  expansion  gave  the  greatest  degree  of  accuracy. 

This  expansion  is 


I(u)  =  a  +  by  +  ( 3 p 2  -1) 


(49) 


Thus,  equations  (46)  and  (47)  yield  the  expression 
3 

l  I(u)  =  a  (  U  +  ”  U_)  +  t^(  M+  2  -  U_2  )  +  ^(u+  3  +  y_3-  p  y_) ( 50) 

k-1 

The  weights  are  found  three  at  a  time  by  solution  of  the  re¬ 
sulting  matrix  equation 
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B .  Material  Properties 


Now  that  the  geometrical  properties  have  been  developed, 
the  physical  characteristics  of  the  radiation  field,  some  of 
which  are  explicitly  angular  dependent,  can  also  be  derived. 

In  this  section,  then,  the  material  properties  of  opacity  and 
emissivity  (in  the  form  of  the  source  function)  and  the  de¬ 
rived  quantity  of  optical  depth  are  developed  in  their  dif¬ 
ference  form. 

Opacity.  Opacity  of  an  extended  atmosphere  is  a  strong 
function  of  radius.  There  are  numerous  methods  used  to  de¬ 
fine  this  dependency  (see  for  an  example  Reference  2).  The 
discretization  formula  chosen  in  this  study  was  a  cubic  spline 
fit,  of  the  form 


X  (r(j))  =  +  a2  r(j)  +  a3  r(j): 


(51) 


where  al ,  a2,  and  a3  are  constants  chosen  to  mimic  the  nature 
of  the  plasma  being  analyzed.  The  validity  of  this  formula 
was  established  by  Dr.  George  Nickel  in  his  cylindrical  ver¬ 
sion  of  the  ray  integration  technique. 

The  Optical  Depth.  From  knowledge  of  the  opacity,  the 
optical  depth  at  the  tangent  points  can  be  derived.  Recall 
that  the  relationship  between  the  two  is 


t  =  JSi  x  (s')  ds' 


Using  the  geometry  established  where 
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(52) 


the  integral  in  ds  can  be  reduced  to  an  integral  in  dr  by 
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The  integral  is  solved  by  the  exact  integration  formula  taken 
from  Reference  1 ,  where 


P  =  /r  2 -r  2 
i  2  o 


r  2-r  2 
i  o 


S3 


P,  =  r 


/r  2 -r  2  -r_  4,  2-r  2  +  i  ( ln(«^r  2-r 


2  2  0  o  1  o 


2  _/r  3_r  2 


2  0  i  o 


P,  =  ( ( r  2  +  2r  2)  /r  2-r  2  -r  2  +  2r  2  /r  2-r  2)/3 

3  2  O  2  °  1  O  i  O 


P4  =  ((r  /r  2  +  2r  2  )  (r  2  +1.5  r  2))/4 

#  2  2  o  2  O 


+  I  ^ro 


(r  -  In  /r  2-r  2) 


-  ((r,  /r,2+  2ro2)  (rt2  + 1.5rQ2))/4 


-  i  <ro  in  •7rrrv>' 


(55) 


The  Source  Function .  The  last  material  property  neces¬ 
sary  for  the  solution  of  the  transfer  equation  is  the  emis- 
sivity.  This  quantity  appears  in  the  transfer  equation  in  the 
form  of  the  source  function,  S. 

It  is  assumed  that  the  source  function,  like  opacity,  is 
a  function  of  radius.  The  exact  variation  is  not  as  simple  to 
predict  as  with  the  opacity.  A  more  complex  relationship  ex¬ 
ists.  Traditionally,  an  estimate  of  the  source  function  is 
made  based  on  observable  quantities.  However,  a  method  of  re¬ 
fining  this  initial  estimate  has  been  produced  and  is  the 
fundamental  quantity  which  makes  the  ray-integration  technique 
both  unique  and  eminently  practical.  In  this  section,  the 
"first  guess"  source  function  estimate  is  derived  and  then, 
from  it,  the  specific  solution. 

Initial  Estimate  of  the  Source  Function.  The  source 
function  of  a  system,  such  as  a  stellar  body  in  the  depths  of 
space  or  a  small  plasma  confined  in  a  magnetic  field  in  a  lab¬ 
oratory,  can  be  inferred  from  the  radiation  flux  emitted  from 
a  spherical  body  of  radius  r  is  just  the  rate  of  energy  change 
per  unit  volume.  Defining  the  quantity  P  (for  power)  as  the 
rate  of  energy  flow,  along  a  distance  r, 

H  =  — l—  (56) 

4tt  r 

where  H  is  the  Eddington  flux  as  defined  in  Chapter  II.  Now 
recall  the  first  moment  equation  is 
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-  XH 


(29) 


3K  +  3K-J 
3r  r 


Assuming  that  the  atmosphere  is  optically  thick  at  r,  then 
according  to  the  Eddington  approximation,  3K  ■  J,  and  the 
second  term  vanishes.  Further,  making  use  of  the  assumption 
of  radiative  equilibrium,  J  =  S,  yields  K  =  1/3  S.  Equation 
(29)  becomes 
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(57) 


Using  a  forward  differencing  scheme  on  the  left  hand  side  and 
assuming  that  the  right  hand  side  can  be  averaged  such  that 


S(i  +  1)-S(.j) 
r(  j  +  1)  -  r(  j) 


=  "I  P 
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Solving  through  for  S(j)  yields 


S(j)  =  S  (  j  +  1 )  +  j  (r(j  +  1)  -  r(j>) 


K(,1)  +  K (  j  +  1 ) 
r  (  j  )  2  r  (  j  +  1 )  2 


(59) 
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Note  that  this  scheme  can  not  be  used  at  the  outer  bound¬ 
ary  since  in  the  steaming  limit  K  =  J,  and  the  second  term 
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does  not  vanish.  Instead,  another  method  must  be  developed 
to  solve  for  the  source  function  at  that  limit.  This  is  done 
by  the  introduction  of  the  extrapolated  boundary.  The  extra¬ 
polated  boundary  is  an  imaginary  boundary  outside  of  the  phys¬ 
ical  limit  of  the  plasma  system  where  the  flux  goes  to  zero. 
The  extrapolated  boundary  is  illustrated  in  Figure  6. 


Figure  6 


Flux  at  the  Boundary  (3:59) 


streaming  limit,  where  by  assumptions  of  Section  II. C. 3  H  =  J. 
Still  assuming  radiative  equlibrium,  J  =  S,  thus  by  equation 
(57), 
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This  is  the  relation  used  to  express  the  source  function  at 
the  approximate  physical  boundry.  The  discrete  form  is 


S(  jn) 


=  JL  X( jn) 
2lT  r(jn) 


(62: 


Last  in  developing  a  source  function  for  the  sphere  is 
the  derivation  of  the  innermost  region.  It  is  assumed  that 
the  centermost  radius  region  contains  a  source  which  increases 
linearly  with  radius  in  the  manner  illustrated  in  Figure  7. 


Within  the  source  region,  the  Eddington  flux  also  varies, 
continuously  rises  from  H(o)  to  some  value  H(2),  as  illus 
trated  by  Figure  8. 


Figure  8 

Radial  Variation  of  the  Eddington  Flux. 


The  Eddington  flux  is  thought  to  vary  linearly  within  this 
region.  The  solution  of  the  first  moment  equation,  thus, 
gives 
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Differencing  the  above  yields 


S(l)  =  S(2 )  +  | 
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Refinement  of  the  Source  Function  Estimate 


the  reasoning  outlined  above,  a  value  for  the  source  can  be 
estimated  at  each  point  on  the  radius  mesh.  However,  recall 
the  solution  to  the  transfer  equation  involved  an  integral  of 
the  source  over  the  length  of  the  ray.  This  equation  is 


I(s2)  =  Ksj )  e 


/s  y (s)ds' 
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To  avoid  having  to  solve  this  rather  cumbersome  integral 
equation,  we  create  an  "equivalent  source  function"  which  we 
will  call  S*,  such  that  equation  (65)  can  be  made  into 


I(s2 )  =  I(Sl  )  e-T  +  S*  (1  -  e  T) 


Development  of  the  S*  Function .  To  establish 
the  S*  function,  first  consider  the  radius  segment  shown  in 
Figure  9. 


Figure  9 


Radial  Segment  Showing  Position  of  r*. 


The  distance  r*  is  defined  as  the  spatial  distance  along  the 


ray  of  length  S0  -  ,  at  which  the  equivalent  source  function 

S*  exists.  Said  in  other  words,  at  r*,  the  source  function 
has  a  value  S*,  which,  if  used  in  the  solution  of  the  transfer 
equation,  reduces  (65)  to  (66).  Assuming  that  the  source 
function  varies  linearly  over  the  interval  S2  to  s^,  the  func¬ 
tion  S*  can  be  approximated  as 


S*  =  S 
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(67) 


Substituting  this  into  (66)  yields 
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Note  that  the  first  integral  term  can  be  easily  reduced  to 


The  second  integral 
as 
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term  can  also  be  reduced  by  defining  r* 
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Using  (69)  and  (70)  in  (66)  gives 
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Thus,  we  have  transformed  the  messy  integral  of  equation  (65) 
into  the  much  neater,  and  numerically  convenient,  form  of 
equation  (66 ) . 


I(s2)  =  I(s,  )  e~T  +  S*  (1-e'T) 


The  difficulty  now  lies  in  the  solution  to  the  r*  inte¬ 
gral.  This  is  done  via  the  method  of  cubic  splines  leading 
to  incomplete  gamma  functions.  The  definition  of  the  method 
of  cubic  spline  interpolation  is  found  in  Appendix  B.  Its 
application  to  generation  of  the  r*  function  is  given  below. 

Development  of  the  r*  Function.  Recall  that 


the  r*  function  is 


fT  r(T ' )  e  T  d  t 


First,  the  interpolatory  function  is  fitted  to  the  r  (t) 


polynomial  which  has  the  values  of 

r  (o)  =  ri 
r  (t)  =  r2 

^lT‘(/r!-ro)/(roX(ro)) 

With  these  values  and  slopes,  r  (t )  is  of  the  form 
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so  that 


Some  care  must  be  taken  when  evaluating  these  incomplete  gamma 
functions  at  small  argument,  since  large  subtraction  errors 
develop.  For  example, 
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As  x->  °°,  r  ->  3!,  as  it  should.  Clearly,  though,  when  t  <<1, 
a  significant  problem  develops  in  that 


r3  (t )  a  £  T1*  (ignoring  e  T) 


If  t  is  10  3,  for  example,  is  2.5x10  ^3.  Clearly,  the  dif¬ 
ference  between  6  and  6(1±10  3  +  1/2x10  1/6x10)  ^e  T  will 

tax  the  decimal  places  of  a  CRAY  in  single  precision. 

However,  the  following  prescription  gives  good  accuracy 
over  the  entire  range. 

Let  Psmall  =  1/4  x^-1/5  t  ^  +  1/12  \ 

Pbig  =  6-(1.02228x3  +  3.00939t2 +  5.002488T +  5.1244)e'  T 
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Then  T3(t)  Psmall  +  Pbig 

The  lower  terms  are  found  by  downward  recursion  such  that 
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(Reference  8)  Thus,  the  r*  integral  is  solved  by  generating 
the  incomplete  gamma  function  in  the  method  described  above  at 
each  ds  interval. 


Numerical  Solution  to  the  Transfer  Equation 

Now,  all  physical  parameters  of  the  plasma  have  been  pre¬ 
sented:  the  dimensions,  the  angular  variables  and  the  material 
properties.  With  this  information  on  hand,  the  intensity  (and 
subsequently,  its  moments),  can  be  solved  for.  Recall,  how¬ 
ever,  that  the  material  properties  input  are  only  approximate. 
Even  the  S*  function  is  based  on  an  approximate  "first  guess" 
source  function,  and  thus  has  limited,  albeit  improved,  accu¬ 
racy.  Consequently,  the  moments  J,  H,  and  K  share  the  limited 
accuracy  of  I.  The  repair  for  this  inaccuracy  is  done  via 
iteration. 

In  this  section,  the  numerical  form  of  the  solution  to 
the  transfer  equation,  that  is,  the  derivation  of  the  value 
for  I  at  each  radius,  is  presented.  Following  this,  the  cal¬ 
culation  of  the  moments  of  intensity  is  expounded  upon.  Last 
is  presented  the  iteration  techniques  used  to  force  the  condi¬ 
tion  of  radiative  equilibrium,  thereby  providing  the  correct 
solution  of  the  problem. 

Solution  of  the  Transfer  Equation .  Using  the  generated 
source  function  S*,  the  transfer  equation  reduces  in  terms  of 
the  radius  mesh  to 


I(j  +  1)  =  I(j)e“T(j  +  1)  +  S*  (1-e  T(j  +  P  (79 


where,  I(J+1)  is  the  intensity  at  the  last  ds  interval.  In 
previous  chapters,  we  have  discussed  passing  rays  through  a 
sphere,  but  have  not  defined  the  order  in  which  this  is  done. 
Figure  10  shows  the  order  of  rays  which  are  evaluated  for  a 
five  radius  radial  mesh. 


Beginning  with  the  outermost  ray,  SI,  the  equation  is  solved 
for  each  segment  along  the  ray  until  the  entire  ray  is  summed 
over,  or  integrated. 

In  this  example,  there  are  five  radii.  The  outer  ray  has 
two  ds  segments  to  the  ray.  The  next  ray  has  four  segments. 
The  last  has  eight.  Beginning  with  ray  one,  1(5)  is  solved 
from 
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where  1(6)  =0.  Then  1(A)  is 


1(4)  =  I(5)e"  (5)  +  S*  (1  -  e'  (5)  ) 


(81) 


Again  1(5)  is  solved  for  by 

1(5)  =  I (4  )e”  (4)  +  S*  (1  -  e'  (4)  )  (82) 
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♦ 

S* 

(1  -  e~ 

(5)> 

—  I (3)  =  I (4  )e” 

(4)  + 

S* 

(1  -  e 

<4)> 

etc . ,  to 

—  1(5)  =  I(4)e~ 

(4)  + 

s* 

(1  -  e~ 

<*>> 

The  process  is  repeated,  ray  by  ray,  until  the  center 
ray  is  integrated. 

Solution  of  the  Moments  of  Intensity .  By  solving  for 
the  intensity  at  each  radius  for  each  ray,  we  are  simply  solv¬ 
ing  for  the  intensity  at  each  angle  on  the  radius  mesh.  The 
moments  of  intensity  are  then  found  by  numerically  integra¬ 
tion,  using  these  solutions.  Or 


J(r)  = 

T 

U  =  ” 1 

I(r,u)  w 

H(r)  = 

T1 

y=- 1 

I  (  r  ,  y  )  w. 

K(r)  = 

M=+» 

I 

I(r»y)  ’ 

U=” 1 
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where  w^  are  the  weights  for  each  angle  and  I  is  the  value  of 
intensity  solved  for  at  each  angle  (ray  interception  point). 

I teration  of  the  Solution.  As  previously  stated,  the  so¬ 
lution  to  the  transfer  equation  has  a  limited  accuracy  as 
solved  above,  due  to  the  initial  crudeness  of  estimates  of 
the  material  properties.  The  source  function  approximates, 
but  does  not  guarantee,  radiative  equilibrium.  However, 
knowing  J  and  S  at  each  radius  point  equilibrium  can  be  forced 
by  successive  iteration.  The  method  of  iteration  depends  im¬ 
plicitly  on  the  optical  thickness  of  the  material.  There  are 
two  commonly  accepted  schemes  of  iteration,  the  Lambda  itera¬ 
tion  scheme  for  the  optically  thin  case  and  the  Unsold  itera¬ 
tion  scheme  for  the  optically  thick  case. 

Lambda  Iteration  Scheme  for  the  Optically  Thin  At¬ 
mosphere  .  The  Lambda  iteration  is  essentially  the  statement 
of  radiative  equilibrium.  That  is,  the  source  at  each  radius 
point  is  simply  replaced  by  the  mean  intensity  calculated  by 
the  ray-integration  scheme.  Analytically, 


S(j)  =  J(j) 


(83) 


I 


It  would  seem  at  first  glance  that  this  is  sufficient  to  es¬ 
tablish  radiative  equilibrium;  however,  it  is  not  over  a  wide 
variety  of  ranges.  As  stated  by  Mihalis  in  Reference  7,  "In 
principle,  successive  applications  of  the  (Lambda  iteration 
scheme)  should  improve  the  solution,  and,  eventually,  produce 
the  exact  solution.  ...In  practice,  however,  the  convergence 


is  too  slow  to  be  of  any  value,  if  the  effective  range  ... 
is  of  order  =  1 ,  so  errors  at  large  depth  are  removed  only 


'infinitely  slowly.'"  (7:62). 

The  Unsold  Iteration  Scheme  for  Optically  Thick 
Atmospheres ■  As  a  result  of  the  Lambda  iteration  being  inap¬ 
propriate  to  use  for  optically  thick  atmospheres,  another 
scheme  must  be  used.  This  second  method  is  called  the  Unsold 
iteration,  after  its  creator  (Ref  7).  The  Unsold  iteration, 
is  developed  from  the  first  order  moment  equation,  which  is 

M  +  =  (27) 

3r  r  * 

Introducing  an  integrating  factor,  q,  known  as  the  sphericity 
function,  the  above  equation  can  be  reduced  to 


v^,v 


1  d  ( q  K ) 
q  dq 


=  "XH 


(84) 


where 


K  d£ 
q  dr 


3K-J 


(85) 


Differencing  this  yields 


tr 


q(j+l)  K(j +1) -q(j)  K(j)« j|x(j)q(j)H(j)+  x(  j+l)q(  j+l)H(  j+1) 


*  (r( j+1) -  r( j)) 


(86) 
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qC  j  +  1  X  (j+l)+7(r(j+i)  -r(j)x(j+l)H(  j+1)^  =  q(j)^x(j)  ~j 

*  (r( j+l)-r( j))x( j)H( (87) 

We  then  assume  that  deviation  from  the  "true"  values  of  the 
moments  vary  linearly  such  that 


1  d(qfkAJ)  =  -  X  (AH) 
q_^_ - 


(88) 


where  fR  =  AK/AJ,  AH  =  (P/A*r  )-Hcalc>AJ  =  Jex  -  Jcalc- 


Note 


that  for  radiative  equalibrium  J  =  S.  We  substitute  the  q's 

calculated  into  the  previous  equation  and  solve  through  for 

J  „(=S  _  ) ,  which  yields 

exact'  exact  ’  7 


S(j)  = 


_  (r(jd)-rQ)K_(j)H(j)  .  (r(.i+l)-r( i) )  P 
J(J)  2ft j7J  +2  ffj)  Wfj)* 

S(.j+l)(PJ(.i+l)/(4irr(  j)2H(  i))+'?'(r(.i+l)  -  r(j))x(i+l)H(.i+l) 
f(j)  +  J(  j+1)  +  j(r(  j+1)  -  r(  j))xCj+l)  +  H  (j+1) 


This  is  a  self  consistent  equation  which,  in  the  limit  of 


AH=0,  produces 


S(j)  =  J(j) 


This  value  can  now  be  used  to  generate  a  new,  more  accurate 
S*  function  and  hence  bring  the  answer  closer  to  the  correct 
equilibrium  value.  This  iteration  solution  brings  about  a 
stable  solution  within  just  a  few  iterations. 


Solution  of  the  Moments  Equations 


■V^ 


>.v 


Although  solution  of  the  moments  equations  will  not 


yield  any  new  information  in  terms  of  the  intensity  through 


the  sphere,  there  is  a  usefulness  in  solving  for  these  equa¬ 


tions  numerically.  That  is,  once  values  for  the  intensity, 


source  function,  J  and  K  have  been  solved  for  they  can  be 


substituted  into  the  moments  equations.  The  moments  equa¬ 


tions  are  now  solved  for  H  and  the  value  of  H  compared  with 


that  generated  by  the  integration  of  rays.  This,  then,  is 


an  excellent  check  on  the  relative  accuracies  of  the  values 


calculated  as  well  as  an  informative  look  at  relative  accu¬ 


racy  of  the  moments  equations  as  compared  to  each  other.  For 


the  purpose  of  completion,  the  dif ferization  of  the  moments 


equations  are  reproduced  below. 


The  Zero' th  Order  Moment  Equation.  The  formulation  for 


the  zero' th  moment  equation  is 


IlllfHi  =  x ( S  -  J) 
r2  3r 


Using  the  backward  difference  method  on  the  term  on  the  left 


hand  side,  and  averaging  the  right  hand  side,  yields 


r( j+1  )2-r( j)2 


r(,j)2H(,j)  r(  j-l)2H(  j-1)  li 


r( j-1)2 


f  X(j)(S(j)-J(j)  + 


X (j+l)(S(j-l)-J( j-1))  1  (92) 


H(0)  is  not  solved  for  in  this  manner ,  but  is  defined  as 
identically  equal  to  zero  for  reasons  outlined  in  Chapter  II 
The  First  Order  Moment  Equation .  In  a  like  manner,  the 
first  order  moment  equation 


£K  +  3K-J 

3  r  r 


is  given  the  numerical  equivalent 


i+1)  -  K(  i)  1  3K(i)-J(i)  3K(.i+l)  -  J(.i+1)  _  1  f  nnr  n 

pTF  r®' 4  1  rTj)  r(ik)  ~  2  x(j)Wj) 

_  J 


x  (j+l)H(j+l)  (94; 


s 


All  of  the  elements  of  the  ray-integration  technique 
have  now  been  discussed.  In  the  first  chapter,  the  basic 
definitions  of  the  physical  nature  of  radiation  transport 
were  examined.  After  which,  an  overview  of  the  technique 
itself  was  provided  along  with  a  development  of  the  logic 
which  spawned  it.  And,  in  the  last  chapter,  the  numerical 
represen ta t ions  of  the  equations  used  in  radiative  transport 
were  presented.  The  efforts  to  so  differentiate  the  equa¬ 
tions  culminate  in  the  creation  of  a  computer  code  which  in 
turn  implements  the  ray-integration  technique  and  produces  an 


in  the  grey  at- 


in  Appendix  D. 


For  convenience  to  the  reader,  the  logic  flow  of  the  program 
is  presented  on  the  following  page. 


Table  VI 


Logic  Flow  for  Computer  Code 


INPUT: 

--  OPACITY 
—  RADII 
--  POWER 


CALCULATE: 

--  SOURCE  FUNCTION 
--  ANGLES  AND  ANGULAR  WEIGHTS 
--  OPTICAL  DEPTH 
--  R*  AND  S* 


FOR  EACH  RAY 

--  FOR  EACH  ds  INTERVAL 

--  SOLVE  TRANSFER  EQ .  FOR  INTENSITY 
--  ACCUMULATE  MOMENTS 


NEXT  RAY 

SOLVE  MOMENTS  EQUATIONS 


OUTPUT 

--  TOTAL  RADIATION  FLUX 
--  NET  RADIATION  FLUX 
--  SOURCE  FUNCTION 
--  EDDINGTON  FACTORS,  f,  and  f 
--  EDDINGTON  FLUX  FROM  K  H 

--  RAY-INTEGRATION  TECHNIQUE 
--  ZERO  AND  FIRST  MOMENT  EQ. 


VI.  Results  and  Discussions 


To  meet  the  second  objective  of  this  thesis,  verifica¬ 
tion  of  the  ray-integration  technique  via  computer  code,  the 
code  stated  in  the  last  section  was  run  with  a  variety  of  in¬ 
puts.  A  sample  of  the  results  produced  and  a  discussion  of 
those  results  are  presented  in  this  section. 

Case  1_ 

To  begin  verification  of  the  program,  the  code  was 
run  with  a  set  of  inputs  approximating  a  sphere  of  outer  ra¬ 
dius  equal  to  0 . Oil  meter ,  constant  opacity,  with  an  input  for 
rate  of  energy  change  per  radius  of  6.28x10^  watts.  The 
radius  mesh  used  consisted  of  22  radii,  evenly  spaced  and  hav¬ 
ing  a  width  of  0.005  meters  between  mesh  points.  These  values 
are  essentially  arbitrary,  in  that  any  physically  valid  set  of 

data  could  be  used.  The  output  to  the  ray-integration  program 

2 

for  this  set  of  data  is  the  total  radiation  flux  (^nr  H),  the 
opacity,  the  source  function,  S,  the  all-angle  radiant  inten¬ 
sity,  J,  and  the  moments  ratios  which  include  the  Eddington 
factors  f^  and  f^  as  well  as  the  ratio  J/S.  The  results  are 
illustrated  in  Figures  11  through  14. 

Consider,  first,  the  total  radiation  flux.  Theory  pre¬ 
dicts  that  the  flux  should  rise  from  zero  in  the  center  to  a 
constant  value  in  the  atmosphere  and,  due  to  the  constraint 
of  radiative  equilibrium,  remain  constant  throughout  the 
sphere.  Figure  11  shows  that  the  program  constructed  follows 
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Opacity  In  unite  of  tnvtm  rtilui 


Opacity 


Source  Function 


Figure  12 

Opacity  and  Source  Function  (Case  I) 


All-angle  Radiant  Intensity 


Eddington  Factors 

1.21  .  .  r.  .  . . 


0.000  0.002 


0.004 


0.006 

Radius 


0.008 


Figure  14 

Net  Radiation  Flux  (Case  I) 
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this  prediction  reasonably  well  in  the  outer  regions  of  the 
sphere,  but  produces  some  fluctuations  in  the  innermost  re¬ 
gions.  Figure  11  shows  the  solution  to  the  flux  for  100  iter¬ 
ations.  The  first  three  iterations  are  identified  on  the 
graph  by  the  numbers  1,  2  and  3.  The  solid  line  which  is 
shown  is  the  value  to  which  the  iteration  converged.  The 
fluctuations  illustrated  are  not  disturbing  since,  in  any  nu¬ 
merical  analysis,  slight  inaccuracies  are  bound  to  occur. 

These  could  be  due  to  round  off  error  by  the  computer  or  to 
non-optimum  approximations  of  the  continuous  functions  in¬ 
volved.  In  this  case,  it  is  believed  that  the  latter  is  re¬ 
sponsible  for  the  fluctuations  which  occur  in  the  output.  To 
support  this  hypothesis,  the  following  explanation  is  given. 

Recall  that  there  are  fewer  rays  intersecting  (and  thus 
fewer  tangent  points  occurring)  at  the  inner  radii  than  at 
the  outer  ones.  Indeed,  only  one  ray  intersects  the  central 
radius.  This  means  that  the  angular  resolution  in  the  center 
regions  is  considerably  less  than  in  the  outer  limits.  Thus, 
any  errors  in  angular  variables  will  be  more  prominent  in  the 
center  regions.  It  is  thought  that  the  errors  are  in  the  an¬ 
gular  weights  approximations.  This  theory  was  examined  by 
varying  the  expansion  of  the  intensity  function  and  genera¬ 
ting  new  angular  weights  for  each  expansion.  The  function 
was  expanded  in  numerous  ways  including  Chebyshev  polynomials 
(Type  I  and  II),  variations  of  expansions  in  terms  of  the 
angle  <p  and  the  sine  of  the  angle  <t>  ,  and  numerous  expansions 


in  terms  of  cosines  of  the  angle  ®.  All  of  these  expansions 


| 

I 


I 

I 
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effected  the  fluctuations  in  flux  by  some  amount:  all  increas¬ 
ing  the  fluctuations.  Only  the  Chebyshev  polynomials,  Type  I, 
were  close  to  the  degree  of  fluctuation  caused  by  the  Lengen- 
dre  expansion,  but  there  was  still  no  improvement.  It  is  be¬ 
lieved  that  there  is  still  an  as  yet  undiscovered  variation  of 
intensity  expansion  which  will  result  in  an  even  better  fit  of 
the  program  output  to  predicted  results. 

Now  consider  some  of  the  other  output  produced  as  shown  in 
Figures  12  -  14.  In  Figure  12,  the  opacity  is  shown  as  a  re¬ 
minder  to  the  programmer  of  the  opacity  input.  This  figure 
also  shows  the  source  function  and  average  intensity  output  at 
each  radius.  The  moments  ratios  (Figure  13)  show  again  how 
reasonably  the  ray-integration  model  obeys  physical  law.  Re¬ 
call  that  for  the  radiative  equilibrium  condition  J=S  in  the 
atomosphere  region,  thus  the  ratio  J/S  should  go  quickly  to 
one.  Recall  also  that  the  Eddington  approximation  states  that 
the  Eddington  factor  f^  should  be  approximately  1/3  at  the  cen¬ 
ter,  (isotropic  limit),  and  progress  toward  1.  In  addition, 
the  ratio  f^  is  expected  to  go  from  0  to  1  over  the  region  of 
the  sphere.  The  results  shown  in  Figure  13  show  good  approxi¬ 
mation  of  these  conditions.  The  quantity  J/S  goes  quickly  to 
one  and  drops  slightly  at  the  end.  This  drop  at  the  end  can 
be  explained  by  recalling  that  an  extrapolated  boundary  was 
imposed  on  the  sphere  which  numerically  created  an  artifically 
high  source  function  at  the  outer  boundary,  thus  resulting  in 
an  artificially  low  J/S  ratio,  which  the  Unsold  iteration  was 
unable  to  correct  for  completely.  Experimentation  with  other 


boundary  conditions  are  in  order  to  rectify  this  condition. 


As  shown  in  Figure  13,  the  ratios  and  begin  at  the 
proper  limits  in  the  center  of  the  sphere,  (1/3  and  0  respec¬ 
tively)  and  build  slowly  toward  one.  The  values  do  not  quite 
reach  one  at  the  physical  limit  of  the  sphere  because  only  be¬ 
yond  this  limit  does  true  streaming  occur  and  the  values  of 
f,  ->  1  and  fu  ->  1  actually  occur. 

Lastly,  consider  the  solution  to  the  moments  equations 
shown  in  Figure  14.  The  solid  line  is  the  value  of  the  Ed¬ 
dington  flux  solved  for  directly  in  the  ray-integration  tech¬ 
nique  and  the  other  lines  were  formed  by  substituting  in  the 
other  moments  calculated  to  the  moments  equation  and  solving 
through  for  H.  Note  that,  in  the  outer  regions,  all  three 
equations  agree  very  well.  In  the  center  regions  of  the  sphere 
the  solution  of  the  zero 1 th  moment  seems  high  and  the  solution 
to  the  first  moment  seems  low.  These  variations  are  also  be¬ 
lieved  due  to  the  inaccuracy  of  angular  weights  in  the  center 
regions  of  the  sphere. 

Case  I I 

The  first  variation  done,  to  experiment  with  the  versa¬ 
tility  of  the  program,  was  to  introduce  a  deliberate  perturba¬ 
tion  to  the  input  source  function.  This  was  done  to  see  if  the 
program  would  still  iterate  to  the  answer  given  in  the  first 
case,  providing  the  mesh,  power  and  opacity  remained  the  same. 
The  perturbation  was  to  introduce  an  artificially  high  central 
source  value  by  multiplying  the  source  created  by  various  val¬ 
ues.  The  total  flux  radiation  output  by  increasing  the  source 
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region  by  a  factor  of  2  appears  in  Figure  15.  The  flux,  which 
the  program  iterated  to,  is  exactly  that  of  the  first  case. 
Several  perturbations  were  introduced,  all  of  which  iterated 
to  the  same  answer.  This  result  indicates  that  the  iteration 
technique  used,  in  all  cases,  brought  about  equilibrium  as  near 
as  it  could,  considering  the  lack  of  angular  resolution. 

Case  III 

Another  variation  done  for  verification  was  to  alter  the 
radius  mesh  of  Case  I,  leaving  all  other  parameters  the  same. 
The  new  radius  mesh  used  also  consisted  of  23  radii,  but  a 
logarithmic  spacing  was  used  in  which  the  spacing  between 
radii  increased  moving  out  from  the  center  to  the  physical 
boundry  of  the  sphere.  The  output  produced  for  the  logarith¬ 
mic  mesh  appear  in  Figures  16  through  19.  Note  that  the  re¬ 
sults  produced  with  the  self-similar  mesh  are  virtually  the 
same  as  for  the  even  mesh  case.  Fluctuations  still  occur 
within  the  inner  regions  of  the  sphere  but,  while  not  as  large 
in  magnitude,  there  are  noticeably  more  of  them. 

This  follows  for  the  error  existing  in  the  angular 
weight  calculated.  That  is,  with  smaller  gaps  between  radii 
in  the  center  of  the  sphere,  the  ds  intervals  integrated  over 
will  be  smaller  and  more  will  exist  closer  to  the  center  of 
the  sphere.  If  the  expansion  used  for  the  intensity  (remem¬ 
ber,  this  generates  the  production  of  angular  weights)  is  not 
accurate,  then  the  departures  of  the  intensity  from  the  true 
value  might  be  more  noticable  over  small  segments  sampled. 

With  the  logarithmic  mesh  used,  there  are  more  samplings  of 
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Figure  17 

Opacity  and  Source  Function  (Case  III) 


AJI-angle  Radian}  Intensity  Eddington  Fodors 


Figure  18 

All-Angle  Radiant  Intensity  and  Eddington  Factors  (Case  III) 


Net  Radiation  Flux 


the  intensity  taken  at  the  inner  radii  of  the  sphere,  thus 
more  fluctuations  will  appear. 

Consider  the  remainder  of  the  output  for  this  case,  Fig¬ 
ures  17  -  19).  First,  note  that  Figure  18  shows  an  improved 
agreement  of  J/S  and  f^.  There  is  still,  however,  fluctua¬ 
tion  in  fH.  Careful  examination  of  the  figure  shows  that  f^ 
does  start  at  zero,  but  rises  quickly  to  approximately  1/4, 
but  then  drops  and  resumes  its  climb  toward  one.  The  sudden 
rise  in  f^  in  the  inner  regions  of  the  sphere  occurs  here 
for  the  same  reason  that  it  did  in  Case  I,  lack  of  angular 
resolution  in  the  center. 

Also  note  that  the  net  radiation  flux  as  illustrated  in 
Figure  19  shows  good  agreement  between  the  different  technique 
used  to  compute  it. 


VII.  Conclusions  and  Recommendations 


£*: 


The  preceding  chapter  has  presented  a  series  of  cases 
for  which  the  ray-integration  technique  was  implemented.  In 
these  cases,  reasonable  agreement  is  made  between  the  results 
produced  and  those  expected  from  the  physics  of  the  situation, 
Inaccuracies  still  exist  in  three  areas.  Most  obviously,  the 
angular  weights  generated  are  not  optimum.  In  addition,  the 
boundary  conditions  used  to  generate  the  source  function  at 
the  outer  boundary  is  not  optimum  in  that  it  produces  an  arti¬ 
ficially  high  value  at  this  limit. 

It  is  suggested  that  further  efforts  be  made  to  find  the 
appropriate  weight  scheme  and  boundary  condition  for  this  pro¬ 
gram  . 

It  is  further  recommended  that  the  technique  developed 
be  expanded  to  applicability  in  a  frequency  dependent  (i.e., 
non-grey  atmosphere)  problem.  Also  expansion  should  be  made 
to  make  the  technique  valid  for  time  varying  energy  flow. 
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Appendix  A 
Gaussian  Quadrature 


Given  a  definite  integral 


j  f(x)  w(x)  dx 
a 


IE  # 


where  f(x)  is  some  polynomial  exisitng  over  the  interval  (a,b) 
This  integral  can  be  approximated  by  the  sum 


n 


l  f(Xk)  Ak 
k=l 

which  contains  2n  +  1  parameters: 

n  xk  1  s,  points  for  evaluating  f(x) 
n  A^  '  s,  coefficients 

and 

the  choice  of  n  itself 

It  was  pointed  out  by  Gauss  that  for  unevenly  spaced  x  points 
for  an  n'th  degree  polynomial,  whichis  orthogonal  to  all  lower 
degree  polynomials,  that  the  A^  coefficients  simply  represented 
a  set  of  weights.  And  that  if  these  weights  were  multiplied  by 
the  function  f(x)  at  the  x  interval,  then  the  summation  approx¬ 
imation  to  the  integral  would  be  exact.  Some  of  the  polyno¬ 
mials  which  fit  these  criterion  are  Legendre,  Laguerre,  Her- 
mite,  and  Chebyshev  (1:970-973). 
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B 


Cubic  Spline  Interpolation 


One  method  of  obtaining  an  estimate  values  for  a  function 
which  can  be  approximated  as  a  polynomial  over  a  large  inter¬ 
val,  is  to  the  divide  interval  into  a  series  of  subintervals 
and  construct  a  polynomial  along  each  subinterval  which  gener¬ 
ally  follows  the  overall  function  in  a  pieceways  manner.  If 
the  polynomial  function  being  interpolated  is  of  third  order, 
then  the  method  used  is  called  cubic  spline  interpolation,  and 
is  described  as  follows. 

Given  a  function  f(x)  defined  on  the  interval  a  to  b 

where  the  interval  is  divided  into  subintervals  x  such  that 

..v  a<x  <  x,  <  . . .  <  x  <  b,  a  cubic  spline  interpolant,  S,  can 

o  i  n 

be  formed  if  the  function,  f,  obeys  the  following  conditions: 

a)  S  is  a  cubic  spline  polynomial,  denoted  S^,  on  the 

subinterval  (x^,  for  each  j  =  0 ,  1,  ...,  n  -  1; 

b)  S(x^)  =  f(xj)>  f°r  each  j  =  0,  1,  ....  n; 

c)  Sj(xj  +  j)  =  sj  +  i^xj  _  for  each  j  =  0,  1,  ..., 

n  -  2 

d)  sj(xj  +  =  sj  +  i (x j  +  i )  for  each  j  =  0,  1,  ..., 

n  -  2 

e)  sj(xj  +  =  s]  +  i(xj  +  x)  for  each  j  =  0,  1,  ..., 

n  -  2 

f)  one  of  the  following  set  of  boundary  conditions  is 
satisfied 

(i)  S"(x  )  =  S" (x  )  =  0 
■„*  o  n 


B-l 


r 


(ii)  S" (xq )  =  f"(x0)  and  S'(xn)  =  f’(xn) 

Using  condition  one  the  interpolant,  S,  is  of  the  form 


S  .  ( x )  =  a.  +  b.h  +  c.h2  +  d.h3 


where  h  =  (x  -  xj)>  if: 


=  f  ( b ) 


=  f  '  ( b ) 

=  3 ( a  -  b)  -  ( f 1 ( a )  +  f'(b)) 
h  h 


2 ( b  -  a)  +  (f ’ (a)  +  f 1 (b)) 
h  h 


(5:116-128) 
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**i 


Term 


al /a2 /a3 
bfac / cfac 
cosphl 
cosph2 
c(  jj) 


delsor 


etaui 


g(  j) 
gs(  j) 


jl/ j2 


hmomO 


hmoml 


nf  req 
phi  (  j  j  ) 
ppr ime 
psml / pbig 
pi  -  4 
sourc ( j ) 


radi 


Appendix  C 


Glossary  of  Computer  Terms 


Definition 


fits  for  opacity  polynomial 

fits  for  cubic  spline  interpolant 

cosine  of  angle  at  beginning  of  ds  interval 

cosine  of  angle  at  end  of  ds  interval 

cosine  of  jj  angle 


opacity  at  j 

slope  of  opacity  at  j 


radius  index 


number  of  radii 
tangent  point  index 

//  of  radius  at  the  beginning/end  of  the  ds 
interval 

solution  for  H  from  zero  moment  equation 
solution  for  H  from  first  moment  equation 
number  of  iterations  to  be  performed 
angle  at  tangent  point  jj 
power /2*pi 

estimates  for  gamma  functions 
exact  solution  of  tau  integrals 
source  function  at  j 
radiant  intensity 


rad j ( j ) 
radh( j ) 
radk(  j ) 


all-angle  radiant  intensity  at  j 
Eddington  flux  at  j 


rstar 


second  moment  of  intensity  at  j 
r*  function 
dr ( jl ) /ds 


IMA 


r  1  /  r2 


dr ( j2 ) /ds 
impact  parameter 

values  of  the  radii  at  the  beginning/end  of  the 
ds  interval 


sstar 


tauO 


tl  -  A 


S*  function 


opacity  along  line  of  sight 
solution  to  tau  integral 


w(  jw) 


weight  at  interval  jw 


APPENDIX  D 


SOURCE  CODE 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 


$OPTIONS  S=600,D=10,A=40,X 

C  - 

C  BEGIN  PROGRAM  BY  INITIALIZING  ARRAYS 

C  - 

DIMENSION  R( 25 ) , WGHT ( 700 ) , TAUE( 10 ) 

DIMENSION  G(25),GS(25), TAU0( 325 ) , SOURC (25) ,C(325 
DIMENSION  PI (325) ,P2(325) ,P3(325) ,P4(325) 
DIMENSION  RADJ ( 25 ) , RADH( 25 ) , RADK( 25 ) , FLUX ( 25 ) 
DIMENSION  DELSOR (25), HMOMO (25), HM0M1 (25) 

REAL  MSQR 

TYPE' THE  PROGRAM  HAS  BEGUN' 

C  - 

C  READ  INPUT  FILE 

C  - 

CALL  OPEN  (5, ’LFILE.TXT", IERROR) 

IF  (IERROR  .NE.  0)  THEN 

TYPE  'CANNOT  OPEN  INPUT  FILE' 
STOP  'RUN  ABORTED' 

ENDIF 

READ  (5,*)  JN,(R(I),I=1,JN),A1,A2,A3,A4,A5,A6 

C  - 

C  PRINT  INPUT  FILE  AND  INITIALIZE  CONSTANTS 

c  - 

TYPE ' JN=  ' , JN 
TYPE ' RADII=  ' 

WRITE  (1,1)  (R(I) ,1=1 , JN) 

1  FORMAT( F7 . 4 ) 

TYPE  ' A1 =  ' , A1 
TYPE  ' A2=  ' ,A2 
TYPE  ' A3=  ' ,A3 
TYPE  ' POWER=  ' ,A4 
TYPE  ' NFREQ=  ',A5 
TYPE  ' START=  ',A6 
START=A6 
P0WER=A4 
NFREQ=A5 
ITER=0 

PI=3. 141592654 
PPRIME=ABS( POWER/ (2*PI)) 

C  - 

C  CALL  SUBROUTINE  TO  CALCULATE  ANGULAR  VARIABLES 


TYPE  'BEGIN  GENERATION  OF  PHI  DEPENDENT  VARIABLES' 
CALL  XFAC  ( JN , R , WGHT ,  C  ,  PI ,P2,P3,P4) 

TYPE  'END  GENERATION  OF  PHI  DEPENDENT  VARIABLES' 

C  - 

C  CALCULATE  OPACITY  AND  SLOPE  OF  OPACITY 


D-l 


51 
5  2 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 


100 

C 

C 

C 


101 

C 

C 

C 

C 

C 

C 

C 

C 


JNM= JN-1 

JNM2=JNM-1 

DO  100  J=1 , JN 

G(J)=(A3*R(J)+A2)*R(J)+A1 

GS ( J ) =2 . 0*A3*R( J ) +A2 


CALCULATE  'FIRST  GUESS' SOURCE  FUNCTION 


SOURC ( JN ) =PPRIME/ (R(JN)*R( JN ) ) 

D0101  JP=1 , JNM2 
J=JN-JP 

SOURC (J )= SOURC (J+l )+(R(J+l)-R(J))*0.75*PPRIME 

*(G(J+1)/(R(J+1)*R(J+1))+G(J)/(R(J)*R(J))) 

CONTINUE 

SOURC ( 1 ) =SOURC ( 2 ) +START*PPRIME*G( 2 ) /R ( 2 ) 


BEGIN  MAJOR  LOOP  TO  CALCULATE  ALL  RAYS 


DO  3000  NT=1 , NFREQ 


CALCULATE  TAU  FOR  EACH  DS  SEGMENT 


DO  300  JZP=1 , JNM 

JZ= JN- JZP 

DO  280  JK=1 , JZP 

JI=JZP*(JZP-l)/2+JK 

J1=JN-JZP+JK-1 

J2=J1+1 

R1=R(J1) 

R2  =  R ( J2 ) 

DR=R2-R1 

DR2=DR*DR 

DR3=DR2*DR 

BB=(3.*G(J2)-3.*G(Jl)/DR2-(GS(J2)+2. *GS ( J1 ) ) /DR 
CC=(GS(J2)+GS(J1))/DR2-2.*(G(J2)-G(J1 ) )/DR3 
T1=(((-CC)*R1+BB)*R1-GS(J1))*R1+G(J1) 
T2=((3.*CC)*R1-2.*BB)*R1+GS(J1) 

T3=-3.*CC*R1+BB 


c 

c 

c 


BEGIN  CALCULATIONS  FOR  ENTIRE  RAY 


m 


1111 


»v 


DO  400  JZP=1 , JNM 
JZ= JN- JZP 
RO=R ( JZ ) 

COSPH2=C( JN*( JN-1 ) /2+JZ 
IF  (NT  .EQ.  100)  WRITE( 1,1111 ) JZ 
FORMAT ( '  JZ= ' 12 / , ' J1  J2  JJ ' , 3X, ' R1 ' ,6X, ' R2 ' , 
7X , ' G1 1 ,10X, ' G2 ' , 6X , 'TAU' ,8X, ' Cl ' ,4X, ' C2 ' , 
7X, ' SI ' ,9X, ' S2 ’ ,9X, ' SSTAR ' ) 

J1 = JN+1 
RADI=0 . 0 


CALCULATE  'INWARD-GOING' 


DO  380  JK=1 ,  JZP 

JJ=JJ+1 

J1=J1-1 

J2=J1-1 

R1=R(J1) 

R2=R( J2 ) 

G1=G( J1 ) 

G2=G(J2) 

JW=J2*(J2-l)/2+JZ 
JI=JZP*(JZP-l)/2+Jl-JZ 
C0SPH1=C0SPH2 
C0SPH2=C( JW) 

TAU=TAU0(JI) 

RS1=C0SPH1/G1 

RS2=C0SPH2/G2 

BFAC=(3.*(R1-R2 ) /TAU- (RSI +2 . *RS2 ) ) /TAU 
CFAC= ( -2.*(R1-R2) /TAU+RS1+RS2 ) /TAU /TAU 
ETAUI=EXP ( -TAU ) 

PSML= ( (0 . 08333333*TAU-0 . 2 ) *TAU+0 . 25 ) *TAU**4 
PBIG=6.-(((1 . 02228*TAU+3 . 00939 ) *TAU+5 . 02488 ) 
*TAU+5 . 1244 )*ETAUI 
TI3=PSML*PBIG/ (PSML+PBIG) 
TERM=TAU*TAU*TAU*TAU*ETAUI 
TI2=0. 333333* (TI3+TERM) 

TERM=TERM/TAU 
TI1=0 . 5*(TI2+TERM) 

TERM=TERM/TAU 
TIO=TI 1 +TERM 

RSTAR=R2+(RS2*TI1+BFAC*TI2+CFAC*TI3)/TI0 
Cl= ( R2-RSTAR )/ (R2-R1 ) 

C2=(RSTAR-R1 )/(R2-Rl) 
SSTAR=C1*S0URC(J1)+C2*S0URC(J2) 

icicieic'k'k’kicie’kie'kicicicic'kicic'kic'kic'kicicic'k'k'k'kie 


TRANSFER  EQUATION 
******************************** 


RADI=RADI*ETAUI+SSTAR*( 1 . -ETAUI ) 


V‘V 


w  v 


153 

154 

155 

156 

157 

158 

159 

160 
161 
162 
163 
!C$ 

165 

166 

167 

168 

169 

170 

171 

172 

173 

174 

175 

176 

177 

178 

179 

180 
181 
182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 

199 

200 
201 
202 

203 

204 


If  (NT  .EQ,  100)  WRITE ( 1,1112)  J1 , J2 , JJ , R1 , R2 ,G1 , 
1  G2 , TAU , Cl , C2 , SOURC ( J1 ) , SOURC( J2 ) , SSTAR 
1112  F0RMAT(2I3,I4,2F7.4,2E11 .2, Ell . 3 , 2F7 . 3 , 3E12 . 4 ) 

C  - 

C  ACUMMULATE  MOMENTS 

c  - 

DRAD=RADI*WGHT( JW ) 

WPERP=-C( JW) 

RADJ ( J2 ) =RADJ ( J2 ) +DRAD 
DRAD=DRAD*WPERP 
RADH( J2)=RADH( J2)+DRAD 
DRAD=DRAD+WPERP 
RADK( J2 )=RADK( J2 )+DRAD 
380  CONTINUE 
JZ1=JZP+1 
JZ2=2*JZP 
J1=JZ-1 

C  - 

C  CALCULATE  'OUTWARD-GOING'  RAY 

C  - 

DO  390  JK= JZ1 , JZ2 

JJ=JJ+1 

J1=J1+1 

J2= J1 +1 

Rl =R ( J1 ) 

R2  =  R( J2 ) 

G1=G( J1 ) 

G2=G( J2 ) 

JW=J2*(J2-l)/2+JZ 
JI=JZP*(JZP-l)/2+J2-JZ 
C0SPH1=C0SPH2 
C0SHP2=C ( JW ) 

TAU=TAUO ( JI ) 

RS1=-C0SPH1/G1 
RS2=-COSPH2 / G2 

BFAC=(3.*(R1-R2)/TAU-(RS1+2.*RS2) )/TAU 
CFAC= ( -2 .*(Rl-R2) /TAU+RS1 +RS2 ) /TAU /TAU 
ETAUI=EXP( -TAU ) 

PSML= ( ( 0 . 083  333  33*TAU-0 . 2 ) *TAU+0 .25) *TAU**4 
PBIG=6.-(  (( 1.02228*TAU+3. 009  39)  *TAU+ 5.02488) 

1  *TAU+5.1244)*ETAUI 

TI3=PSML*PBIG/ (PSML+PBIG) 

TERM=TAU*TAU*TAU*ETAUI 
TI2+0. 333333*(TI3+TERM) 

TERM=TERM/TAU 
TI1=0 . 5*(TI2+TERM) 

TERM= TERM /TAU 
TI0=TI1 +TERM 

RSTAR=R2+(RS2*TI1+BFAC*TI2+CFAC*TI3)/TI0 
Cl = ( R2-RSTAR )/ (R2-R1 ) 

C2  =  (RSTAR-R1)/ (R2-R1  ) 

SSTAR=C1 * SOURC (Jl ) +C2*S0URC ( J2 ) 


D-4 


TRANSFER  EQUATION 
RADI =R ADI *ETAUI +SSTAR * ( 1 . -ETAUI ) 

'kic'k'k’k-k'kic'k'k'k'kic'kicic'k'k'k'k'k'k’k'k'k'klc'k'k'k'k'k 

IF  (NT  .EQ.  100)  WRITE  (1,1112)  J1 , J2 , JJ , R1 , R2 , 
G1 , G2 , TAU , Cl , C2 , SOURC ( J1 ), SOURC (J2 ), SSTAR 

ACUMMULATE  MOMENTS 


DRAD=RADI*WGHT( JW) 

WPERP=C ( JW ) 

RADJ ( J2 ) =RADJ ( J2 ) +DRAD 
DRAD=DRAD*WPERP 
RADH( J2 ) =RADH( J2 ) +DRAD 
DRAD=DRAD*WPERP 
RADK( J2 ) =RADK( J2 ) +DRAD 
CONTINUE 

GO  BACK  AND  CALCULATE  NEXT  RAY 

RADH( 1 )=0 . 0 
RADK( 1 ) =RADK( 1 ) / 3 . 0 
CONTINUE 

DELSOR( 1 )= SOURC ( 1 )-RADJ ( 1 ) 

DO  440  J=2 , JN 

DELSOR ( J )  =  SOURC ( J ) -RADJ ( J ) 

DEFINE  AND  PRINT  FLUX  FOR  ALL  RADII 

FLUX(J)=4*PI*(R(J)*R(J))*RADH(J) 
CONTINUE 
RADH( 1 ) =0 . 0 
FLUX( 1 ) =0 . 0 

IF  (NT  .EQ.  100)  GO  TO  441 
ITER=ITER+1 
TYPE 

TYPE  'ITERATION  '.ITER 
TYPE  'TOTAL  RADIATION  FLUX  RADIUS' 
DO  557  1=1, JN 

WRITE  (1,557)  FLUX( I ) , R( I ) 

FORMAT  ( El 1 . 3 , 30X , F7 . 4 ) 

CONTINUE 

IF  (POWER  .GT.  0.0)  GO  TO  445 

DO  LAMBDA  ITERATION  FOR  THIN  CASES 

SOURC ( 1 )=0 . 0 
DO  443  J=2 , JN 
SOURC (J )=RADJ ( J) 

GO  TO  449 
CONTINUE 


mmfm 


257:  C  - 

258:  C  DO  UNSOLD  ITERATION  FOR  THICK  CASES 

259:  C  - 

260:  SOURC ( JN )  =  SOURC ( JN ) *  P PR I ME/ (2.*R(JN)*R(JN) 

261:  1  /RADH(JN) 

262:  DO  446  JP=1,JNM2 

263:  J=JN-JP 

264:  C0N=(R(J+1 )-R(J) )*0. 5 

265:  A=CON*G(J)*RADJ(J)/RADK(J) 

266:  B=CON*G( J  +1 ) *RADJ ( J  +1 ) / RADK ( J  +1 ) 

267:  SOURC (J )=(RADJ( J )-A*RADH( J ) )*(SOURC( J+l )+B*PPRIME/ (2. 

268:  1  R(J+l)*R(J+l)))/( RADJ (J+l ) +B*RADH( J+l ) ) +A*P PRIME/ 

269:  2  (2.*R(J)*R(J)) 

270:  446  CONTINUE 

271:  B=A 

272:  CON=R ( 2 ) *0 . 5 

273:  A=CON*G ( 1 )  *RADJ ( 1 ) /RADK ( 1 ) 

274 :  SOURC ( 1 ) = ( RADJ ( 1 ) -A*RADH( 1 ) ) * ( SOURC ( 2 ) +B*PPRIME/ ( 2 . * 

275:  1  R(2)*R(2)))/( RADJ ( 2 ) +B*RADH( 2 ) ) 

276:  SOURC ( 1 ) =SOURC ( 1 ) + ( PPRIME/ ( 2 . *R( 2 ) *R ( 2 ) *RADH( 2 ) ) ) * 

277:  1  ( SOURC ( 1 ) -RADJ ( 1 ) ) 

278:  449  CONTINUE 

279:  3000  CONTINUE 

280:  C  - 

281:  C  SOLVE  FIRST  AND  ZERO ' TH  MOMENT  EQUATIONS 

282:  C  - 

283:  HMOMO ( 1 ) =0 . 0 

284:  HMOMl ( 1 ) =0 . 0 

285:  DO  455  J=2,JNM 

286:  RBAR=0.5*(R(J-1)*R(J-1)+R(J)*R(J)) 

287:  RDIFF=R(J)-R(J-1) 

288 :  RHSBAR=0. 5*(G( J-l )*DELSOR( J-l )+G( J)*DELSOR( J) ) 

289:  HMOMO (J)=R(J-1)*R(J-1) *HMOMO (J-1)/(R(J)*R(J)) 

290:  1  +RDIFF*RBAR*RHSBAR/(R(J)*R( J)) 

291:  455  HMOMl (J )=-(( RADK( J+l )-RADL( J-l ))/( R( J+l ) -R( J-l ) ) 

292:  1  + ( 3 . 0*RADK( J ) -RADJ (J))/R(J))/G(J) 

293:  HMOMl (JN)=-( ( RADK( JN)-RADK( JN-1 ) ) / ( R( JN) -R ( JN-1 ) ) 

294:  1  + ( 3 . 0*RADK( JN )-RADJ(JN) )/R(JN) )/G(JN) 

295:  C  - 

296:  C  PRINT  OUTPUT 

297:  C  - 

298:  TYPE 

299:  TYPE 

300:  TYPE  'OPACITY  RADIUS' 

301:  DO  149  1=1, JN 

302:  WRITE  (1,150)  G(I),R(I) 

303:  150  FORMAT  ( Ell . 3 , 7X , Ell . 3 ) 

304:  149  CONTINUE 
305:  TYPE 

306:  TYPE 

307:  TYPE  'SOURCE  RADIUS' 

308:  TYPE' FUNCTION' 

309:  DO  151  1  =  1, JN 


310 

311 

312 

313 

314 

315 

316 

317 

318 

319 

320 

321 

322 

323 

324 

325 

326 

327 

328 

329 

330 

331 

332 

333 

334 

335 

336 

337 

338 

339 

340 

341 

342 

343 

344 

345 

346 

347 

348 

349 

350 

351 

352 
253 

354 

355 

356 

357 

358 

359 

360 

361 

362 


WRITE  (1,150)  SOURC(I) ,R(I) 

151  CONTINUE 
TYPE 
TYPE 

TYPE  'ALL  ANGLE  RADIUS' 

TYPE  'RADIANT  INTENSITY' 

DO  152  1=1, JN 

WRITE  (1,150)  RADJ ( I ) , R ( I ) 

152  CONTINUE 
TYPE 
TYPE 

TYPE  ' *************EDDINGTON  FACTOR*************** 
TYPE 

TYPE'  J/S  H/J  K/J  RADIUS' 

DO  500  J=1 , JN 
R ADH ( J ) =RADH ( J ) / RADJ ( J ) 

R ADK ( J)= RADK ( J ) / RADJ ( J ) 

TFAC=0 . 0 

IF  (SOURC(J)  .NE.  0.0)  TFAC=1 . /SOURC( J ) 

RADJ ( J ) =TFAC*RADJ ( J ) 

WRITE  (1,499)  RADJ ( J ) , RADH ( J ) , RADK ( J ) , R ( J ) 

499  FORMAT  ( 3(E11 . 3 ,4X) ,F7 . 4 ) 

500  CONTINUE 
C 

C 

TYPE 

TYPE 

TYPE' *********  NET  RADIATION  FLUX  ************' 
TYPE 

TYPE 'RAY  INTEG  ZERO  MOM  FIRST  MOM  RADIUS' 

DO  505  J  =  1  , JN 

RADH( J )=RADH( J )*SOURC( J ) 

WRITE  (1,251)  RADH ( J ) , HMOMO ( J ) , HM0M1 ( J ) , R ( J ) 

251  FORMAT  ( 3(E11 . 3 ,4X) ,F7 . 4) 

505  CONTINUE 
END 

c  - 

C  BEGIN  GEOMETRY  FACTOR  SUBROUTINE 

C  - 

$OPTIONS  S=200 , X 

SUBROUTINE  XFAC ( JN , R , WGHT , C , PI , P2 , P3 , P4 ) 

DIMENSION  PHI (32 5) ,W(325) ,R(25) , WGHT (700) ,C(700) 
DIMENSION  P1(250),P2(250),P3(250),P4(250) 

REAL  MSQR 
PI-3.14159265 


CALCULATE  PHI  ANGLES  AND  ANGULAR  WEIGHTS 


C(l)=1.0 
C ( 2 ) -1 . 0 
C( 3 ) =0 . 0 
C ( 4 )  =  1 . 0 


C 

C 

C 


363 

36A 

365 

366 

367 

368 

369 

370 

371 

372 

373 

374 

375 

376 

377 

378 

379 

380 

381 

382 

383 

384 

385 

386 

387 

388 

389 

390 

391 

392 

393 

394 

395 

396 

397 

398 

399 

400 

401 

402 

403 

404 

405 

406 

407 

408 

409 

410 

411 

412 

413 

414 

415 


C ( 5 )=SQRT (l.-((R(2)*R(2))/(R(3)*R(3)))) 

C ( 6 )=0 . 0 
WGHT( 1 ) =1 . 0 
WGHT ( 2 ) =1 . / 6. 

WGHT( 3 ) =2 . / 3 . 

WGHT(4)=-0. 099955 
WGHT(5)=0. 3556495677 
WGHT(6)=0. 4889186828 
J  J  =  7 

JNM= JN-1 
DO  1000  KK=4 , JN 
KMAX=2*KK-1 
DO  4  K=1  ,  KMAX 
4  W ( K ) =0 . 0 

KKP=KK+1 
DO  10  K=1 , KK 
TESTX=R(K)/R(KK) 

TESTY=SQRT( 1 . -TESTX*TESTX ) 

IF  (TESTY  .NE.  0.0)  GO  TO  23 
PHI (K) =90.0* .0174532925 
GO  TO  24 

23  ARG=TESTX/ TESTY 
PHI (K)=ATAN(ARG) 

24  CONTINUE 
10  CONTINUE 

DO  20  K=KKP , KMAX 
20  PHI (K)=PI-PHI (KMAX-K+1 ) 

DO  100  K=2 , KMAX , 2 
FM=PHI ( K-l ) 

FO=PHI ( K ) 

FP=PHT(K+1) 

CP=COS ( FP ) 

CO=COS ( FO ) 

CM=COS ( FM) 

DF=CM-CP 

DS=(CM*CM-CP*CP)/2. 

DS2=(CM*CM*CM-CP*CP*CP-CM+CP)/2. 

PSQR=(3.*CP*C:-1. )/2. 

0SQR=(3.*C0*C0-1. )/2. 

MSQR=(3.*CM*CM-1. )/2. 

DEN= ( CO*PSQR-CP*OSQR ) + ( CP*MSQR-CM*PSQR) 

1  +(CM*OSQR-CO*MSQR) 

DWN= ( CO*PSQR-CP*OSQR ) *DF+ ( OSQR-PSQR )*DS+(CP-C0)*DS2 
DWO= ( C  P*MSQR-CM*PSQR ) *DF+ ( PSQR-MSQR ) *DS  +  ( CM-CP ) *DS2 
DWP= ( CM*OSQR-CO*MSQR ) *DF+ ( MSQR-OSQR ) *DS+ ( CO -CM ) *DS2 
W(K-1)=W(K-1 )+DWM/DEN 
W( K ) =DWO /DEN 
W  ( K  +  l )=W(K+1 )+DWP/ DEN 
100  CONTINUE 

DO  110  K=1 ,KK 
C(JJ)=COS(PHI(K) ) 

WGHT(JJ)=W(K)/2. 

no  jj=jj+i 


A 1 6  :  1000  CONTINUE 

417:  TYPE' FINISH  GENERATION  OF  COSINES  AND  WEIGHTS' 

418:  C  - 

419:  C  BEGIN  GENERATION  OF  P1-P4  CONSTANTS 
420:  C 

421:  TYPE' BEGIN  GENERATION  OF  P1-P4  CONSTANTS' 

422:  DO  2000  JZP=1,JNM 

423:  JZ= JN- JZP 

424:  JKl=JZP*(JZP-l)/2+l 

425:  JK2=JK1+JZP-1 

426:  DLAST=0 . 0 

427:  ELAST=0 . 0 

428:  IF  (JZ  .NE.  1)  ELAST=ALOG( R ( JZ ) ) 

429:  DO  2000  JK=JK1,JK2 

430:  J1=JZ+JK-JK1 

431:  J2=J1+1 

432:  Rl=R( J1 ) 

433:  R2=R( J2 ) 

434:  RO=R( JZ ) 

435:  DNOW=SQRT (R2*R2-RO*RO) 

436:  DD=R2+DN0W 

437:  ENOW=0 . 0 

438:  IF  (DD  .GT.  0.0)  ENOW= ALOG ( DD ) 

439:  PI ( JK ) =DN0W-DLAST 

440:  P2 ( JK )  =  ( R2  *DN0W-R1 *DLAST+R0*R0*R0* ( ENOW-ELAST ) ) *0 . 5 

441 :  P3(JK)=((R2*R2+2.0*R0*R0 ) *DN0W- (R1*R1 +2 . 0*R0*R0 ) 

442:  1  *DLAST(*0. 33333 

443:  P4 ( JK )=0 . 25*R2*DNOW*(R2*R2+l . 5*R0*R0)+0 . 375 

444:  1  *RO*RO*RO*RO*ENOW-0 . 25*R1*DLAST 

445:  2  *(R1*R1+1 . 5*R0*R0)-0. 375*RO*RO*RO*ELAST 

446:  DLAST=DNOW 

447:  ELAST=ENOW 

448:  2000  CONTINUE 

449:  RETURN 

450:  END 
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